NDEigensystem and radial function equation for Hydrogen atom

Question

I'm trying to numerically solve the radial equation for the 3D hydrogen atom problem, i.e., to find $R(r)$ which satisfies: $$ -\frac{\hbar^2}{2m}\left[\frac{1}{r}\frac{d}{dr}\left(r^2\frac{dR(r)}{dr}\right)-\frac{l(l+1)}{r^2}R(r)\right]-\frac{1}{4\pi\epsilon_0 r}R(r)=ER(r). $$

The problem is that the NDEigensystem gives me non-sense answers. This is my code:

1. First, I set all the constants as the unity: $\hbar=m=\epsilon_0=1$ and $l=0$.

h = 1;
m = 1;
ϵ0 = 1;
Z = 1;
e = 1;
l = 0;
a0 = (4 π ϵ0 h^2)/(m e^2);

2. I define the Hamiltonian:

Hcoul = 
  -(h^2/(2 m))*(D[R[r], {r, 2}]/2 + 
    1/r D[R[r], r] - (l (l + 1))/r^2 R[r]) - (Z*e^2)/(4 π ϵ0 r) R[r];

3. I use the PDEigensystem routine as follows:

{vals, funs} = 
  NDEigensystem[
    {Hcoul, DirichletCondition[R[r] == 0, True]}, R[r], {r, 0, 2000}, 10, 
    Method -> 
      {"Eigensystem" -> {"Arnoldi", "Criteria" -> "RealPart"}, 
    "SpatialDiscretization" -> 
      {"FiniteElement", {"MeshOptions" -> {"MaxCellMeasure" -> 0.05}}}}];

From the above, I get the following eigenvalues:

{2.89232*10^-6, 0.0000188806, 0.0000364341, 0.0000554983, 
0.0000760327, 0.0000980063, 0.000121395, 0.000146177, 0.000172338, 
0.000199864}

but analytically I know that the answer is $$ E_n=-\frac{1}{32\pi^2n^2}=\left\{-0.00316629, -0.000791572, -0.00035181, -0.000197893, -0.000126651, \ -0.0000879524\right\}. $$

Moreover, when I plot the numerical wave functions, comparing with the analitycal solution:

f[r_, n_, l_] := 
  Sqrt[((2 Z)/(n a0))^3*(n - l - 1)!/(2 n ((n + l)!)^3)]
    Exp[-((Z r)/(n a0))] ((2 Z r)/(n a0))^l 
      LaguerreL[n - l - 1, 2 l + 1, (2 Z r)/(n a0)];

Show[
  Plot[{f[r, 1, 0]}, {r, 0, 200}, 
    PlotStyle -> {{Dashed, Blue}}, PlotRange -> All], 
  Plot[Evaluate[funs[[1]]], {r, 0, 500}, 
    PlotRange -> All, PlotStyle -> Blue]]

I get this:

The dashed curve is the analytical solution whereas the line is numerical. Additionally, many numerical solutions are the same:

Show[Plot[Evaluate[funs], {r, 0, 500}, PlotRange -> All]]

Do you know what I'm doing wrong? Is there something that I'm not considering?

Answer 1

First, we can put $\frac{e^2}{4 \pi \epsilon _0}=1, \hbar =1, m=1 $ and define Hamiltonian in the standard form as

H = 1/2 (-(D[R[r], {r, 2}] + 2/r D[R[r], r]) + l (l + 1)/r^2 R[r] - 
     2/r R[r]);

For this Hamiltonian we know exact solution in a form $$E_n=-\frac {1}{2 n^2}, n=1,2,...\\ R_{ln}=c\rho ^le^{-\rho/2}L^{2l+1}_{n+1}(\rho ), \rho=\frac{2r}{n}$$ Second, we can't get the right solution with numerical method like FEM since for this method we need to limit $r$ while in the real problem Hamiltonian H defined on infinite interval. But we can get some eigenvalues using next code

{vals, funs} = 
  NDEigensystem[{H}, R[r], {r, 0, 100}, 20, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", 
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}}];

En = Take[Sort[vals], 7]

Out[]= {-0.125, -0.0555556, -0.03125, -0.02, -0.0139179, -0.0107609, -0.00732718}

We can compare En to the analytical solution

Table[-1./(2 n^2), {n, 2, 8}]

Out[]= {-0.125, -0.0555556, -0.03125, -0.02, -0.0138889, -0.0102041, -0.0078125} 

We also can compare funs and $R$ for given En as follows

f[r_, n_, l_] := 
  Exp[-((r)/(n ))] ((2 r)/(n ))^l LaguerreL[n - l - 1, 
    2 l + 1, (2 r)/(n )];
ff = Take[
  SortBy[Table[{vals[[i]], funs[[i]]}, {i, Length[vals]}], First], 7]

Table[Plot[{Abs[1.5 f[r, n, 0]/f[0.0001, n, 0]], 
   Abs[ff[[n - 1, 2]]]}, {r, 0, 10}, PlotRange -> All, 
  PlotLabel -> -.5/n^2], {n, 2, 8}]

Therefore funs and $R$ are correlated well with sufficient normalization on absolute value.

Answer 2

Thanks, @AlexTrounev.

I had worked on the problem and I found a mistake in my code: the Laplacian operator was wrong.

Now, your code works very well, but it does not reproduce de ground state energy $E_{n=1}=-0.5$.

With my old strategy and by correcting the laplacian, both codes give the same results and I can get the ground state. This is my current program:

1. Define the constants and Hamiltonian:

h = 1;
m = 1;
ϵ0 = 1;
Z = 1;
e = 1;
l = 0;
a0 = (4 π ϵ0 h^2)/(m e^2);

H1 = -(h^2/(
    2 m))*(D[R[r], {r, 2}] + 
     2/r D[R[r], r] - (l (l + 1))/r^2 R[r]) - (Z*e^2)/( 
   4 π ϵ0 r) R[r] + s*R[r]

I had added a factor $s R(r)$ inspired by the answer given here. It helps to ensure all the eigenvalues are positive definite.

2. I proceed to compute the eigenvalues:

s=10;
{valsH3, funsH3} = 
  NDEigensystem[{H1, DirichletCondition[R[r] == 0, True]}, 
   R[r], {r, 0, 100}, 10, 
   Method -> {"Eigensystem" -> {"Arnoldi", "Criteria" -> "RealPart"}, 
     "SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" -> \
{"MaxCellMeasure" -> 0.001}}}}];

Finally, I obtain:

In[]:= valsH3 - s

Out[]= {-0.499009, -0.124876, -0.0555188, -0.0312345, -0.019992, -0.0138638, -0.00958567, -0.00464891, 0.00166775, 0.00926097}

Now by using your code and the extra factor $sR$ the ground state energy is also obtained:

{valsH1, funsH1} = 
  NDEigensystem[{H1}, R[r], {r, 0, 100}, 20, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}}];

In[]:= Take[Sort[valsH1] - σ, 7]

Out[]= {-0.5, -0.125, -0.0555556, -0.03125, -0.02, -0.0139179, \
-0.0107609}

The wave functions are related to a normalization constant (blue is the exact solution and red the numerical one):

F[r_, n_, l_] := 
  Sqrt[((2 Z)/(n a0))^3*(n - l - 1)!/(2 n ((n + l)!)^3)]
    Exp[-((Z r)/(n a0))] ((2 Z r)/(n a0))^
   l LaguerreL[n - l - 1, 2 l + 1, (2 Z r)/(n a0)];

GraphicsGrid[{Table[
   Plot[Abs[F[r, i, 0]], {r, 0, 10}, PlotRange -> All], {i, 1, 5}], 
  Table[Plot[Abs[Evaluate[funsH3[[i]]]], {r, 0, 10}, PlotRange -> All,
     PlotStyle -> Red], {i, 1, 5}]}, Spacings -> 15, Frame -> All]