How to increase the integration domain of NDEigensystem without a "non-Hermitian" error?

Question

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.

Thank you for your help.

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 https://mathematica.stackexchange.com/questions/198455/using-multiple-boundary-conditions-with-ndeigensystem/198507#198507 — 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 05:19
@AlexTrounev, what is wrong with the usage of NDEigensystem here? — user21, May 20 '19 at 05: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 What is there in your equation, that there may appear a chaotic spectrum? Only numerical instability is there. — Alex Trounev, May 20 '19 at 11:26
@AlexTrounev, OP is not specifying inhomogeneous DirichletConditon or NeumannValue. The note you refer to is to say that only homogeneous DirichletCondition or NeuamnnValue can be used. It does not say that you can not use PeriodicBoundaryCondition. In fact there are examples in the ref page that show the usage of PeriodicBoundaryConditon with NDEigensystem — user21, May 20 '19 at 12:37
@AlexTrounev, I have clarified that in the note section by explicitly stating that PeriodicBoundaryCondition can be used. — user21, May 20 '19 at 12:45
@user21 Of course, we can use PeriodicBoundaryCondition. Moreover, the OP wants to get numerical chaos. — Alex Trounev, May 20 '19 at 14:51
@AlexTrounev, Sorry for the late reply. I suppose it isn't really clear that one expects a chaotic eigenspectrum from this differential equation. It's related to Arithmetic Quantum Chaos (a good paper: http://web.math.princeton.edu/sarnak/Arithmetic%20Quantum%20Chaos.pdf). I expect the spacing of the eigenvalues to be Poissonian in nature. — sr101studios, May 21 '19 at 19:43
@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

user21 · Answer 1 · 2019-05-20T05:43:26.987

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["NDSolve`FEM`"]
(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 05:20

How to increase the integration domain of NDEigensystem without a "non-Hermitian" error?

1 Answers1