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.

$\endgroup$
  • $\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$ – Alex Trounev May 18 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$ – sr101studios May 19 at 5:19
  • $\begingroup$ @AlexTrounev, what is wrong with the usage of NDEigensystem here? $\endgroup$ – user21 May 20 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$ – Alex Trounev May 20 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$ – Alex Trounev May 21 at 21:11
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.

$\endgroup$
  • $\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$ – sr101studios May 21 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$ – user21 May 22 at 5:20

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.