1

As advised in the previous two questions (Why, if I enter an angle into the function, then does the code not work correctly?, let us solve the problem for the Coulomb potential in spherical coordinates

From the general theory, we know the levels:

The ground state — 1s (m=0) = -0.5
The first excited level is fourfold degenerate — 2s (m=0), 2p (m=-1), 2p (m=0), 2p (m=1) = -0.125
The second excited level is ninefold degenerate — 3s (m=0), 3p (m=-1), 3p (m=0), 3p (m=1), 3d (m=-2), 3d (m=-1), 3d (m=0), 3d (m=1), 3d (m=2) = -0,056
...

ClearAll["Global`*"]
rmax = 20;

VS[r_] := -1/r
{valss, funss} = 
  NDEigensystem[{(-1/2*
        Laplacian[ψ[r, θ, ϕ], {r, θ, ϕ}, 
         "Spherical"] + VS[r]ψ[r, θ, ϕ]) + ψ[
       r, θ, ϕ]0.5}, ψ[r, θ, ϕ], {r, 0, 
    rmax}, {θ, 0, Pi}, {ϕ, 0, 2*Pi}, 100, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", 

{"MeshOptions" -> {"MaxCellMeasure" -> 0.05}}}, 
     "Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}];
Sort[valss] - 0.5


{-0.5035, -0.22311, -0.125494, -0.125039, -0.125014, 

-0.081024, -0.0804358, -0.0803342, -0.0617875, -0.060486, -0.0604803, 

-0.0580498, -0.0580485, -0.0580056, -0.0480634, -0.0480223, 

-0.0479788, -0.0461506, -0.0461338, -0.0461152, -0.0371368, 

-0.0371171, -0.0371139, -0.0370261, -0.0333038, -0.0333024, 

-0.0332551, -0.0286157, -0.0286054, -0.0285738, -0.0284433, 

-0.0269414, -0.0269313, -0.0227607, -0.0199048, -0.0198853, 

-0.0198145, -0.0198087, -0.0195684, -0.0154964, -0.0154651, 

-0.0154304, -0.0107441, -0.0107237, -0.0106342, -0.0105755, 

-0.0102346, -0.00118218, -0.00109607, -0.00102089, -0.00100278, 

-0.0009178, -0.000733773, -0.000719368, -0.000136159, 0.00266871, 

0.0027101, 0.00271696, 0.0029014, 0.00931309, 0.00940575, 0.0094768, 

0.00975244, 0.00984893, 0.0102233, 0.0105106, 0.0203588, 0.0205096, 

0.0205552, 0.0206402, 0.020662, 0.0207215, 0.0207282, 0.0210018, 

0.0210577, 0.0212097, 0.0224455, 0.0251264, 0.0251286, 0.0252013, 

0.0321901, 0.03224, 0.0325208, 0.032655, 0.0332945, 0.0334383, 

0.0344648, 0.0386394, 0.0386791, 0.0388232, 0.038835, 0.0393229, 

0.044794, 0.0448802, 0.0452931, 0.0453073, 0.0453284, 0.0453843, 

0.0454258, 0.0457998}

The first eigenvalue agrees with theory (within accuracy) = -0.5035. The next eigenvalue is -0.22311 but this level should not be. Next, we see a threefold degenerate level with the correct eigenvalues -0.125494, -0.125039, -0.125014, but it should be degenerate fourfold. Further, again, the eigenvalues -0.081024, -0.0804358, -0.0803342, -0.0617875, -0.060486, -0.0604803 do not correspond to the theory. The next three eigenvalues roughly correspond to the theoretical values, but again there are three -0.0580498, -0.0580485, -0.0580056, instead of nine.

Why doesn't the code return the correct eigenvalues?
How to rewrite the code so that Mathematica returns the correct levels of eigenvalues (taking into account degeneration)?

Mam Mam
  • 1,843
  • 2
  • 9
  • Yes, this is the way to tackle the problem. At least we see that we get degenerate states as is expected from general theory which you did not get previously because you used cylindrical coordinates. As a matter of fact you should get EXACTLY the numbers exspected from the analytic results. So the exact formulation needs to be checked for errors. One thing is to increase rmax and see if the results get more precise. – Michael Weyrauch Jul 10 '23 at 12:24
  • 1
    Also there might be an issue arrising due to the fact that the Arnoldi Eigenvalue solver orders the eigenvalues according to absolute size (So signs may not be correct). This can be fixed by shiftiting the spectrum and then remove the shift. You have a shift of .5: why? – Michael Weyrauch Jul 10 '23 at 12:24
  • 1
    Maybe you also need to change the discretization mesh to get better results. – Michael Weyrauch Jul 10 '23 at 13:04
  • You should consult also this thread to get more insight: https://mathematica.stackexchange.com/questions/105298/how-to-numerically-solve-a-1-d-time-independent-schr%C3%B6dinger-equation – Michael Weyrauch Jul 10 '23 at 13:12
  • @Michael Weyrauch, thanks for the comments! I do not quite agree with you that degenerate states are obtained. I would say that some additional states have turned out that are superfluous and do not correspond to the theory. I have checked different radii and meshs, it doesn't help, anyway Mathematica gives wrong extra values. – Mam Mam Jul 11 '23 at 08:36
  • 1
    Did you see that in the thread I linked to above the user @Jens gets exactly the correct energy values using NDEigensystem? And what do you mean with extra states??: a atomic p state (in Schrödinger theory) is exactly 3 times degenerate. A straighforward diagonalization method should find this degeneracy. – Michael Weyrauch Jul 11 '23 at 20:58

0 Answers0