# How to increase the integration domain of NDEigensystem without a “non-Hermitian” error?

Recently, I've been using NDEigensystem to solve a two-dimensional eigenfunction equation. The region that I'm solving over is as follows.

keyhole = ImplicitRegion[x^2 + y^2 >= 1, {{x, -1/2, 1/2}, {y, 0.0001, ymax}}];


It's a rather specific region and my boundary conditions are very particular to this region but what seems to be an issue is tuning the ymax parameter. When I set ymax to 100, for example, the NDEigensystem function works perfectly and gives me great results. However, when I increase ymax over 100, I keep getting the following error about being non-Hermitian:

Eigensystem::herm: The matrix SparseArray[Automatic,<<2>>,{1,{{0,<<955>>},{<<1>>}},{0.0167943,0.000175333,0.0000158879,-2.49834*10^-11,5.12298*10^-9,-2.49834*10^-11,0.0000158879,-1.26762*10^-9,1.09711*10^-8,-1.28059*10^-9,4.37737*10^-13,-7.94393*10^-6,4.37737*10^-13,1.09711*10^-8,<<755244>>,-0.000465693,0.0000205407,0.0000581381,0.000116049,-6.88058*10^-19,0.0000291058,-0.000193762,-3.68362*10^-18,4.28843*10^-23,-1.3667*10^-6,-2.92675*10^-23,-0.0000102905,0.0232857}}] is not Hermitian or real and symmetric.


I may be making a silly syntax error with the code but I'm really unsure. The code I'm using is as follows.

test = NDEigensystem[{-y^2 (D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}]),
PeriodicBoundaryCondition[u[x, y], x == -1/2,
TranslationTransform[{1, 0}]],
PeriodicBoundaryCondition[u[x, y], x^2 + y^2 == 1,
Function[x, {{-1, 0}, {0, 1}}.x]]},
u[x, y], {x, y} \[Element] keyhole, 30];


I'd really appreciate your input and feedback about what I can do to increase ymax without invoking this error.

• This is because you do not want to use NDEigensystem[] correctly. If you use the code that I recommended to you, then there are no problems with the growth of ymax see mathematica.stackexchange.com/questions/198455/… – Alex Trounev May 18 '19 at 21:57
• @AlexTrounev, Sorry for not responding to your earlier comment. It's just that the condition that you recommended (that u[x,y]==0 on the boundaries) does not give me the eigenspectrum that I'm trying to compute. Essentially, if I compute it with more relaxed periodic boundary conditions, the spectrum should be chaotic (Poissonian) because I'm trying to compute the eigenspectrum of a particle travelling on a special manifold. I'm very sorry if I wasn't very clear about that in the previous post. I would love to hear your comments on how to resolve this issue. – sr101studios May 19 '19 at 5:19
• @AlexTrounev, what is wrong with the usage of NDEigensystem here? – user21 May 20 '19 at 5:32
• @user21 Read tutorial: Homogeneous DirichletCondition or NeumannValue boundary conditions may be included. When no boundary condition is specified on part of the boundary \[PartialD]\[CapitalOmega], then this is equivalent to specifying a Neumann 0 condition. – Alex Trounev May 20 '19 at 11:10
• @sr101studios Thanks for the preprint. Why did you decide that we need periodic boundary conditions for the realization of the desired spectrum? – Alex Trounev May 21 '19 at 21:11

You'd need to play a bit with different meshes. Here is a refined mesh that works:

ymax = 110;
keyhole =
ImplicitRegion[x^2 + y^2 >= 1, {{x, -1/2, 1/2}, {y, 0, ymax}}];
test = NDEigensystem[{-y^2 (D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}]),
PeriodicBoundaryCondition[u[x, y], x == -1/2,
TranslationTransform[{1, 0}]],
PeriodicBoundaryCondition[u[x, y], x^2 + y^2 == 1,
Function[x, {{-1, 0}, {0, 1}}.x]]},
u[x, y], {x, y} \[Element] keyhole, 30,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.125}}}];


Also note that I changed the y starting point to be 0.

Probably a better way to do to is to create a mesh that is refined at the curve. Use a ymax=10 to see the effect. Part of the issue, I assume, is that the mesh is long and stretched:

ymax = 110;
keyhole =
RegionDifference[Rectangle[{-1/2, 0}, {1/2, ymax}], Disk[]];
Needs["NDSolveFEM"]
(mesh = ToElementMesh[keyhole, AccuracyGoal -> 5])["Wireframe"]
test = NDEigensystem[{-y^2 (D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}]),
PeriodicBoundaryCondition[u[x, y], x == -1/2,
TranslationTransform[{1, 0}]],
PeriodicBoundaryCondition[u[x, y], x^2 + y^2 == 1,
Function[x, {{-1, 0}, {0, 1}}.x]]},
u[x, y], {x, y} \[Element] mesh, 30];


If you want to push it further you could try to make a triangle mesh around the curvature and a quad mesh in the upper part and merge the two meshes.

• Thank you so much for your help! Your code definitely resolves the issue and allows me to increase the region of integration. I have one more question about normalisation. At the moment, the default normalisation computes integral( psipsi )dx dy and normalizes accordingly. However, I want to change the default normalization to integral (psi psi)dx dy / y^2 = 1 as my normalisation condition. Is there a way of doing this using the Method->{"VectorNormalization"} condition? Thank you for your help. – sr101studios May 21 '19 at 19:51
• @sr101studios, yes, that should be possible but that is a new question. Start with something simple and experiment with different normalizations. – user21 May 22 '19 at 5:20