4
$\begingroup$

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);
  1. 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];
  1. 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

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]]

Plot of all numerical eigenfunctions

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

$\endgroup$
5
$\begingroup$

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}]

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

$\endgroup$
3
$\begingroup$

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.

  1. 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]

Comparison between exact solution (blue) and numerical (red)

$\endgroup$
3
  • $\begingroup$ If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review $\endgroup$
    – MarcoB
    Dec 1 '20 at 22:49
  • 1
    $\begingroup$ I wrote the last message as a complement to the answer given by @AlexTrounev. Should I edit it? $\endgroup$ Dec 1 '20 at 23:00
  • 1
    $\begingroup$ @JorgeCastaño Thank you very much for your answer. Normally it could be just update to you main post. It is very good you got the ground state. Actually I know this trick with adding $sR$ and have used it on this forum as well. $\endgroup$ Dec 1 '20 at 23:24

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.