# NDEigensystem cannot solve numerically the 3D Coulomb problem, while DSolve returns the right answer

After having derived by hand the eigenvalues and eigenfunctions for the 3D and 2D hydrogen atom, I want to solve the systems numerically using Mathematica. I need to do this because my next step is to solve numerically some Coulomb-like systems with no analytical solution, and so I need to convince myself that Mathematica is interpreting the problem correctly and giving me reasonable calculations.

The differential equation I am trying to have Mathematica solve is the radial equation for the Coulomb potential in 3D: $$R'' + \frac{2}{r}R' + 2\left[E + \frac{1}{r} - \frac{l(l+1)}{2 r^2} \right]R = 0$$

The solution (radial eigenfunction) for this ODE is of the form:

$$R_{n,l}(r) = \left( \frac{r}{n} \right) ^{l} e^{-\frac{r}{2 n}} L_{n-l-1}^{2l+1} \left( \frac{r}{n} \right)$$ where $$L_n^k(r)$$ is the generalized Laguerre polynomial.

Mathematica knows how to solve this equation symbolically with no boundary conditions specified and on a half-infinite space, $r \in [0,\infty)$:

DSolve[{y''[x] + (2/x)*y'[x] + (-(1/4) - ((l*(l + 1))/(x^2)) + (n/x))*
y[x] == 0}, y[x], x]


yields $$\left\{\left\{y(x)\to c_1 e^{l \log (x)-\frac{x}{2}} U(l-n+1,2 l+2,x)+c_2 e^{l \log (x)-\frac{x}{2}} L_{-l+n-1}^{2 l+1}(x)\right\}\right\}$$ And then after the fact I know to disregard the first term in the solution because of the boundary condition that the function must not diverge at infinity.

Now, when I try to solve the same differential equation using NDEigensystem, I get results that are very wrong. First let me state that I am aware that the argument of NDEigensystem only takes the Hamiltonian operator from the Schrodinger equation, or specifically that NDEigensystem takes only the LHS of $\hat{H}\psi = E\psi$, while for DSolve I am writing the ODE in the form $\left( \hat{H} - E \right) \psi = 0$.

I'll screenshot my code and results below, for both $l=0$ and $l=1$.

{vals3DL0, funs3DL0} = NDEigensystem[-(R''[r]/2) - (R'[r]/r) - (R[r]/r), R[r], {r, 0, 200}, 2]
Plot[Evaluate[funs3DL0], {r, 0, 200}]


And for $l=1$,

{vals3DL1, funs3DL1} = NDEigensystem[{-(R''[r]/2) - (R'[r]/r) + (2*R[r]/(r^2)) - (R[r]/r), DirichletCondition[R[r] == 0, True]}, R[r], {r, 0, 200}, 2]
Plot[Evaluate[funs3DL1], {r, 0, 200}]


EDIT for clarity: the closed form solutions for the eigenfunctions I am expecting are given above as $R_{n,l}$. The eigenfunctions are implemented and plotted in Mathematica below.

Reig[r_, n_, l_] := ((r/n)^l) * Exp[-r/(2*n)]*LaguerreL[n - l - 1, 2*l + 1, r/n];
Plot[Evaluate[Table[Reig[r, n, 0], {n, 1, 4}]], {r, 0, 80}, PlotRange -> {{0, 80}, {-1, 2}}]


As you can see, the solution I am getting from Mathematica oscillates wildly and does not appear to tend exponentially to zero.

• Could you clarify a bit what eigenvalues and/or functions you would have expected? And possibly remove the graphics and put mathematica code there? – user21 Apr 25 '16 at 0:56
• I think you should plot the correct eigen functions so that it's more clear. – xslittlegrass Apr 25 '16 at 1:28
• It's worth pointing out that the l=0 case was the topic of a previous question. You could probably adapt my answer using the coordinate transformation to this case. – Jens Apr 25 '16 at 1:59
• I specified the eigenfunctions that I expected when i wrote the expression $R_{n,l}(r)$. I will add a plot of these eigenfunctions per xslittlegrass' request. Jens, I'll take a look at your previous response (there's a lot to digest) and see if it works. – Matthew Brunetti Apr 25 '16 at 6:35

Your main problem is that you're getting somewhat correct results, but fail to properly interpret them. Help for NDEigensystem reads:

NDEigensystem[<...>] gives the n smallest magnitude eigenvalues and eigenfunctions for the linear differential operator ℒ over the region Ω.

Note the smallest magnitude here. But if you're looking for algebraically smallest eigenvalues, i.e. the closest ones to the ground state, then you must shift the Hamiltonian's spectrum to make ground state the smallest in magnitude. You can do this by adding a constant ($-1$ would suffice, since ground state $E=-n^{-1}=-1$) to your expected eigenvalue $E$ (or, equivalently, negative of that constant to the potential). In NDEigensystem notation this corresponds to addition of R[x] to the first argument.

Then you'll get something much closer to the expected solution:

{vals3DL0, funs3DL0} = NDEigensystem[-R''[r] - (2 R'[r])/r - (2 R[r])/r + 1 R[r], R[r],
{r, 0, 200}, 2]
Plot[Evaluate[funs3DL0], {r, 0, 200}, PlotRange -> All]


{{-0.0408090089390994, 0.798242974843909}, {InterpolatingFunction[{{0., 200.}}, {...

Now, this is not exactly nice result. But this is just because of coarse grid used for the calculation. If you follow the help page on the Options of NDEigensystem, you'll get something like this:

{vals3DL0, funs3DL0} = NDEigensystem[-R''[r] - (2 R'[r])/r - (2 R[r])/r + R[r],
R[r], {r, 0, 200}, 2, Method -> {"SpatialDiscretization" ->
{"FiniteElement", {"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}}]
Plot[Evaluate[funs3DL0], {r, 0, 200}, PlotRange -> All]


{{-1.98185450817448*10^-7, 0.7499999750288},...

I am not sure I understand the question 100% but here is what I think you are looking for:

{vals3DL0, funs3DL0} =
NDEigensystem[-(R''[r]/2) - (R'[r]/r) - (R[r]/r), R[r], {r, 0, 200},
3, Method -> {"Eigensystem" -> {"Arnoldi",
"Criteria" -> "RealPart"},
"SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" -> \
{"MaxCellMeasure" -> 0.05}}}}];


This then gives:

vals3DL0
(* {6.33602*10^-6, 0.00252879, 0.00536387} *)


Eigenvalue (and thus NDEigenvalue) can have a bit of a odd way to sort the result. You could sort according to the RealPart. Also, you may want more elements if one edge is at infinity.

• By "more elements" do you mean I should extend my domain to > 200? This is the dimensionless hydrogen atom so I am effectively already simulating out to 200 Bohr radii. Changing the upper bound on $r$ does not change my eigenfunctions at all. Regarding the eigenvalues that you calculated, I read those results as saying the ground state has the eigenenergy $E_1 = 6.33 * 10^{-6}$, the first excited state has energy $E_2 = 0.00252..$, etc. This must be incorrect because the ground state energy for the hydrogen atom is large, and excited states are closer to zero (more weakly bound). – Matthew Brunetti Apr 25 '16 at 6:42