3
$\begingroup$

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.

Improve this question
$\endgroup$
10
  • $\begingroup$ 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/… $\endgroup$ Commented May 18, 2019 at 21:57
  • 1
    $\begingroup$ @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. $\endgroup$ Commented May 19, 2019 at 5:19
  • $\begingroup$ @AlexTrounev, what is wrong with the usage of NDEigensystem here? $\endgroup$ Commented May 20, 2019 at 5:32
  • $\begingroup$ @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. $\endgroup$ Commented May 20, 2019 at 11:10
  • 1
    $\begingroup$ @sr101studios Thanks for the preprint. Why did you decide that we need periodic boundary conditions for the realization of the desired spectrum? $\endgroup$ Commented May 21, 2019 at 21:11

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Commented May 21, 2019 at 19:51
  • $\begingroup$ @sr101studios, yes, that should be possible but that is a new question. Start with something simple and experiment with different normalizations. $\endgroup$ Commented May 22, 2019 at 5:20

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.