How to numerically solve a 1-d time-independent Schrödinger equation?

Question

The point is to solve the eigensystem of the given Hamiltonian.

I tried ParametricNDSolve combined with FindRoot to search for eigenvalues as xslittlegrass did in his post Find eigen energies of time-independent Schrödinger equation in the first place. It gave almost perfect results,but then I found it very inconvenient for large amounts of eigenstate as xslittlegrass did.

Apparently some of the methods used in xslittlegrass's post won't work here directly. My potential is not a polynomial. I'm trying to adjust their methods but so far no luck.

I tried different approaches, including using the radial wavefunction of hydrogen atom as basis and change the Hamiltonian into a matrix, using finite difference method to convert it into a eigenvalue problem and so on. I then found out the new version of Mathematica added a new function called NDEigensystem.

I tried those all and the outcome was just far away from what I expected. The closest one was the radial wavefunction of hydrogen atom, and it had about 10% error. The others didn't give me a proper result, not even close. As for NDEigensystem, I don't know how to apply a Neumann condition to it. For harmonic oscillator, NDEigensystem seems fine, but when I applied it to Coulomb potential, the results were terribly wrong. I supposed it's because of the improper boundary, I only attached a DirichletCondition to it. But I couldn't add a Neumann condition.

Here's some details (with Jens' help, I've been able to apply NDEigensystem to this, and I'm very grateful for that):

My potential is:

V[r_] = -(1/r) - (1.0415223038416566` E^(-0.9990999998788636` r))/r

1. Using the radial function of hydrogen atom, I wrote the following code:

   a = 1; l = 0; nmax = 20;
   R[n_, r_] = FullSimplify[
      2/(a^(3/2) n^2 (2 l + 1)!) Sqrt[(n + l)!/(n - l - 1)!]
       E^(-1/2 (2 r)/(n a)) ((2 r)/(n a))^l 
         Hypergeometric1F1[-n + l + 1, 2 l + 2, (2 r)/(n a)]]
   u[n_, r_] = r R[n, r]
   HMax = ParallelTable[NIntegrate[u[n, r] (-α^-1 \!\(
     \*SubscriptBox[\(∂\), \(r, r\)]\((u[i, r])\)\) + 
      V[r] u[i, r]), {r, 0, ∞}, 
      AccuracyGoal -> ∞, WorkingPrecision -> 30, 
      PrecisionGoal -> 6], {n, 1, nmax}, {i, 1, nmax}];
   N[Eigenvalues[HMax], 6]
   

   And the results:

   $$ \small\begin{array}{lllll}\{&-1.00113,& -0.15272,& -0.0619985, &-0.0336547, &-0.0211352,\\ &-0.0145071,& -0.0105749,& -0.00805096,& -0.00633443, &-0.00511407,\\ &-0.00421542, &-0.00353452,& -0.00300626,& -0.00258818, &-0.00225162,\\ &-0.00197665, &-0.0017491, &-0.00155863, &-0.00139753,& -0.00125987\quad\}\end{array}$$

   The results,as mentioned above,were not so good.

2. This is the code using NDEigensystem as Jens showed (I'm only interested in eigenvalues,so I replaced NDEigensystem with NDEigenvalues, and I modified it a little to increase accuracy):

   AbsoluteTiming[
    With[{shift = 10, d = 500, n = 20}, 
     ev = NDEigenvalues[{shift f[r] + V[r] f[r] - 
        1/2 f''[r](*+NeumannValue[0,r\[Equal]d]*), 
        DirichletCondition[f[r] == 0, True]}, f[r], {r, 0, d}, n, 
        Method -> {"SpatialDiscretization" -> {"FiniteElement", \
        {"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}, 
        "Eigensystem" -> {"Arnoldi", MaxIterations -> Infinity}}];
     evShifted = ev - shift]]
   

   It took me 2540.17s to complete the calculation.

   $$\small\begin{array}{lllll}\{&-1.33732, &-0.186434, &-0.0710575,& -0.0374072, &-0.023051, \\&-0.0156184, &-0.0112778, &-0.00852433, &-0.00666881, &-0.00535929,\\ &-0.0044008, &-0.00367822, &-0.00312003, &-0.00267982,& -0.0023228,\\ &-0.00199518, &-0.00162513, &-0.0011884, &-0.000689036,& -0.000132268\quad\}\end{array}$$

   As I mentioned in the comments,in low energy levels, the results are good,but when it comes to higher levels, the accuracy drops rapidly. The accuracy only increased with d and MaxCellMeasure, and the code becomes less efficient in the meantime.

   Oh, BTW, a MaxCellMeasure up to 0.001 could be unnecessary. When I change it to 0.01, and increase b to 1000, I spent 1511.33s and have a much better result. For error tolerance only at 10E-7, I suppose bigger MaxCellMeasure would be better.

And here's the binding energy I obtained from the NDSolve & FindRoot approach

$$\small\begin{array}{cccc} n & \text{energy} & n & \text{energy} \\ 1S & 1.33732 & 6S & 0.0156184 \\ 2S & 0.186434 & \dots & \\ 3S & 0.0710575 & 10S & 0.00535929 \\ 4S & 0.0374072 & 20S & 0.0012937 \\ 5S & 0.023051\end{array}$$

It's very kind of Jens to modify this table for me, which was so crude in the first place. The eigenvalue should be like that, but in low energy states, my results are far from what it should be. I don't know where this error is coming from, and if there is a better way to do this.

If this problem can be fixed, I'd be very grateful for that.

Solving Coulomb Potential using Jens' variable substitution method

I'm going to save myself some trouble, just post my code and result here:

code which is basically copy from Jens':

V[r_] = -14.3996448504451555789912633575295935414`6.88249108029756/r;
α = 0.26246843532124502683678435500618410715`6.076375074285658;

radialEq = V[r] f[r] - f''[r]/α;
radialξ = 
   Simplify[radialEq /. f -> (ψ[ArcTan[#]] &) /. r -> (Tan[ξ]),
      Pi/2 > ξ > 0]

With[{max = 20, shift = 10}, {ev, ef} = 
 NDEigensystem[{radialξ + shift ψ[ξ], 
   DirichletCondition[ψ[ξ] == 0, 
   True]}, ψ[ξ], {ξ, 0, Pi/2}, max, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", \
   {"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}, 
   "Eigensystem" -> {"Arnoldi", MaxIterations -> 40000}}];
 evNew = ev - shift]

And results:

$$\small\begin{array}{lllll} &-13.6057,& -3.40142,& -1.51174, &-0.850356,& -0.544228, \\&-0.377936, & -0.277667, &-0.212588, &-0.167969, &-0.13605,\\& -0.112427, &-0.0944449, \ & -0.0804237, &-0.0692623, &-0.0601369, \\&-0.0524875, &-0.0468007, \ & -0.0402104, &-0.0350137, &-0.0310644\end{array}$$

And the real energy:

Ehydrogen = 
 Table[-QuantityMagnitude[UnitConvert[Quantity[1, "Rydbergs"], "eV"]]/
   n^2, {n, 20}]

$$\small\begin{array}{lllll}&-13.60569,& -3.401423, &-1.511744, &-0.850356, &-0.5442277, \\&-0.3779359, \ &-0.2776672, &-0.2125889,& -0.1679715, &-0.1360569,\\& -0.1124437, \ &-0.0944840, &-0.0805071, &-0.06941680,& -0.06046974, \\&-0.05314724, \ &-0.04707852, &-0.04199288, &-0.03768890, &-0.03401423\end{array}$$

The error:

$$\small\begin{array}{lllll}&-2.72446*10^{-8}, &-6.8044*10^{-9}, &-2.97381*10^{-9}, &-7.99033*10^{-10}, & 6.76493*10^{-9},\\& 4.52332*10^{-8}, &2.05556*10^{-7}, &7.58432*10^{-7}, & 2.39862*10^{-6}, &6.70562*10^{-6}, \\&0.0000169127, &0.0000390625, &0.0000833608, \ & 0.000154477,& 0.000332808,\\& 0.000659734,& 0.00027787,& 0.00178247, \ & 0.00267516,& 0.00294987\end{array}$$

When I increase the MaxCellMeasure to 0.0001, the error decrease a lot:

$$\small\begin{array}{lllll}&-2.75312*10^{-8}, &-6.81183*10^{-9}, &-3.03523*10^{-9}, &-1.70471*10^{-9}, \ &-1.08962*10^{-9}, \\&-7.52408*10^{-10}, &-5.35404*10^{-10}, &-3.49364*10^{-10}, \ &-9.34864*10^{-11},& 4.12964*10^{-10}, \\&1.53102*10^{-9}, &3.96303*10^{-9}, \ &9.01172*10^{-9},& 1.89823*10^{-8}, &3.77832*10^{-8}, \\&7.17968*10^{-8}, \ &1.31103*10^{-7},& 2.31154*10^{-7}, &3.9495*10^{-7},& 6.58026*10^{-7}\end{array}$$

Answer 1

Why the original matrix approach fails

The question originally showed an attempt at a solution based on converting the differential operator (the Hamiltonian) into a matrix (HMax) by forming a Table of overlap integrals. The functions used in these integrals were the bound-state eigenfunctions of the hydrogen radial equation. Although the matrix obtained in this way has eigenvalues that converge quickly, the convergence is not to the expected eigenvalues of the potential V[r] (which differs from that of the hydrogen atom).

The reason why this approach doesn't work is that the bound states of hydrogen alone do not form a complete set in which an arbitrary radial wave function can be expanded. What's missing is the uncountably infinite set of positive-energy solutions to the hydrogen radial equation. Their form is in fact given on Mathematica's help page for Hypergeometric1F1, under "Applications." These functions have to be included when considering matrix elements of an arbitrary radial operator. However, since they form a continuum, they don't lend themselves to the construction of a simple matrix as was intended in the question.

The conclusion is that we cannot expect the original approach in the question to work on physical grounds. The numerical approach based on finite-element discretization discussed below circumvents this problem to a certain extent. However, the presence of a continuum in the spectrum indirectly causes trouble in the NDEigensystem calculations, too, because the bound states we're after crowd together in a dense sequence of eigenvalues near the ionization threshold at zero eigenvalue, if you seek progressively larger number of radial nodes.

To salvage the above method of constructing a finite matrix, one should choose a different basis in which to expand the Hamiltonian. For example, the harmonic oscillator in 3D gives you a complete set of radial functions with a purely discrete spectrum that could be used instead of the hydrogenic solutions. Of course, this then becomes more of a physics problem, not a Mathematica issue.

Using NDEigensystem

The question mentions that you didn't know how to use NDEigensystem because of the boundary conditions. So I'll try to address the calculation with this approach. You can in fact use NDEigensystem for this problem, as long as the states of interest fit well into the simulation domain. Here is something that works for the lowest states. It imposes DirichletCondition at both ends. The reduced radial wave function should indeed vanish for $r=0$, too. This is because it's divided by r when forming the 3D wave function, and there can't be a divergence in that function.

Below, I'll choose $\alpha=2$ because it gives agreement with the eigenvalues you were expecting (except for the sign, which I'm pretty sure is simply omitted in the question):

V[r_] = -(1/r) - (1.0415223038416566` E^(-0.9990999998788636` r))/r

With[{shift = 10, d = 30, n = 5}, {ev, ef} = 
  NDEigensystem[{shift f[r] + V[r] f[r] - 1/2 f''[r],
     DirichletCondition[f[r] == 0, True]}, f[r], {r, 0, d}, n, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}, 
     "Eigensystem" -> {"Arnoldi"}}];
 evShifted = ev - shift]

{-1.33732, -0.186434, -0.0710474, -0.0348373, -0.00119517}

With[{n = 4, d = 30, amplitudes = {-1, 5, 20, 30}},
 Plot[Evaluate[
   Table[evShifted[[i]] + amplitudes[[i]] 1/r ef[[i]], {i, n}]], {r, 
   0, d}, PlotRange -> {-4, 4}, 
  Epilog -> {Gray, Dashed, 
    Table[Line[{{0, evShifted[[i]]}, {d, evShifted[[i]]}}], {i, n}]}]]

Here I plotted the wave functions of the first four states (divided by r because that's how your f[r] is related to the true 3D radial wave function). The functions are offset by their eigenvalues. Looking at the list of eigenvalues above, you can also see that the agreement with the values listed in the question is good, for the first three states.

The fourth wave function with the largest number of nodes is seen to extend all the way to the boundary where I applied the DirichletCondition. This causes numerical deviations from the true eigenvalue of the infinite domain for all states higher than the third one. You can fix this by increasing the domain and changing MaxCellMeasure.

In order to get the eigenvalues to come out sorted the right way, I manually added a shift to the potential such that all eigenvalues are positive. Then I subtract that shift again at the end. It's supposed to be possible to specify the same method for the eigensolver with the Arnoldi method, but I found that adding the method option "Arnoldi ", "Shift" -> -10 while setting my manual shift to zero, the eigenvalues are sorted the wrong way around for our purposes of getting the spectrum starting with the ground state (so it's mainly a matter of taste whether to use that option or not).

Edit more efficient calculation

Using NDEigensystem on the given radial equation is made difficult by the fact that the integration domain extends from zero to positive infinity. Since the solver is based on a finite-element discretization, I decided to reformulate the problem using a transformation of variables that maps the radial coordinate r onto a finite interval: $r\equiv \tan\xi$ where $\xi \in [0,\frac{\pi}{2}]$. This produces accurate results much faster than the more straightforward approach above.

First, Let's do the substitution of variables in the radial equation (assuming V is defined as above):

Clear[f,r,ξ,ψ];
radialEq = V[r] f[r] - f''[r]/2;

radialξ = 
 Simplify[radialEq /. f -> (ψ[ArcTan[#]] &) /. 
   r -> (Tan[ξ]), Pi/2 > ξ > 0]

$$-\frac{1}{2} \cos ^4(\xi ) \psi ''(\xi )+\sin (\xi ) \cos ^3(\xi ) \psi '(\xi )+\psi (\xi ) \left(-1.04152 e^{-0.9991 \tan (\xi )}-1\right) \cot (\xi )$$

This is now the equation to be solved, and the rest proceeds as before (again I use only Dirichlet boundary conditions at both ends):

With[{max = 20, shift = 10}, {ev, ef} = 
  NDEigensystem[{radialξ + shift ψ[ξ], 
    DirichletCondition[ψ[ξ] == 0, True]}, ψ[ξ], {ξ, 0, Pi/2}, max, 
    Method -> {"SpatialDiscretization" -> {"FiniteElement", 
     {"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}, 
      "Eigensystem" -> {"Arnoldi", MaxIterations -> 40000}}];
 evNew = ev - shift]

$$\begin{array}{cccc} -1.33732 & -0.186434 & -0.0710575 & -0.0374072 \\ -0.023051 & -0.0156184 & -0.0112777 & -0.0085241 \\ -0.00666806 & -0.00535715 & -0.00439514 & -0.00366612 \\ -0.00310643 & -0.00260488 & -0.00223437 & -0.00192196 \\ -0.00162614 & -0.0014017 & -0.000978152 & -0.000545353 \\ \end{array}$$

The agreement with the last results listed in the question (using NDEigenvalues with large domain) seems to be very good, but the speed is orders of magnitude faster with the transformation of variables as done here.

Edit: comparison to hydrogenic wave functions

The hydrogen atom is a good test case because it exhibits the same divergence at the origin as the potential in the question, and also is very long-range so that the trick of transforming the independent variable is equally useful here. The exact solutions are known, and I define them below as u for the special case of zero angular momentum (that's what the radial equation in the question refers to).

Then I repeat the numerical solution for the Coulomb potential and plot a comparison of the solutions with the exact results in a single panel for all solutions:

Clear[u]; 
u[n_][r_] := 
 r (2 E^(-(r/n)) Sqrt[(-1 + n)!/n!] LaguerreL[-1 + n, 1, (2 r)/n])/n^2

Clear[f, r, ξ, ψ];
V[r_] := -1/r;
radialEq = V[r] f[r] - f''[r]/2;
radialξ = 
 Simplify[radialEq /. f -> (ψ[ArcTan[#]] &) /. r -> (Tan[ξ]),
   Pi/2 > ξ > 0]

$$ -\frac{1}{2} \cos ^4(\xi ) \psi ''(\xi )+\sin (\xi ) \cos ^3(\xi ) \psi '(\xi )-\psi (\xi ) \cot (\xi )$$

This is the equation to be solved next. Note that I ask for the solutions to be returned as bare InterpolationFunction objects by omitting the argument of ψ[ξ] in the second slot of NDEigensystem. This makes the subsequent manipulations easier. Other than that, I used the same settings as in the previous calculation.

With[{max = 20, shift = 10}, {ev, ef} = 
   NDEigensystem[{radialξ + shift ψ[ξ], 
     DirichletCondition[ψ[ξ] == 0, True]}, ψ, {ξ, 0,
      Pi/2}, max, 
    Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}, 
      "Eigensystem" -> {"Arnoldi", MaxIterations -> 40000}}];
  evNew = ev - shift];

Needs["NumericalCalculus`"]

amplitudes = 
  Table[NLimit[u[i][r]/ef[[i]][ArcTan[r]], r -> 0], {i, Length[ef]}];

GraphicsGrid@
 Partition[
  Table[Plot[{amplitudes[[i]] ef[[i]][ArcTan[r]], u[i][r]}, {r, 0, 
     30}, Frame -> None, Axes -> None, PlotRange -> All, 
    PlotStyle -> {Red, Directive[Black, Dashing[.1]]}], {i, 
    Length[ef]}], 4]

The two plots lie on top of each other for all 20 states.

To make the comparison, I had to find the proper amplitude factors that make the graphs equal at some point. I chose to do that by calculating the ratio between numerical and exact solution in the limit $r\to 0$ using NLimit. This limit always exists even though both functions vanish there.

The eigenvalues cluster together around zero for higher node numbers, and therefore a constant absolute error translates into a much larger relative error. This can only be fixed by increasing the spatial resolution. For the (low) resolution that I chose in the above calculation, the eigenvalues are only useful up to roughly the tenth solution:

evNew[[;; 14]]

$$\begin{array}{cccc} -0.5 & -0.125 & -0.0555556 & -0.03125 \\ -0.02 & -0.0138889 & -0.010204 & -0.00781215 \\ -0.00617175 & -0.00499706 & -0.00412515 & -0.0034518 \\ -0.00293943 & -0.00250752 & \text{} & \text{} \\ \end{array}$$

-0.5/Range[14]^2

$$\begin{array}{cccc} -0.5 & -0.125 & -0.0555556 & -0.03125 \\ -0.02 & -0.0138889 & -0.0102041 & -0.0078125 \\ -0.00617284 & -0.005 & -0.00413223 & -0.00347222 \\ -0.00295858 & -0.00255102 & \text{} & \text{} \\ \end{array}$$

The last line shows the expected exact results.