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

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.39964485044515557899126335752959354146.88249108029756/r;
α = 0.262468435321245026836784355006184107156.076375074285658;

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

With[{max = 20, shift = 10}, {ev, ef} =
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}$$

• You might find this useful: Solving the Schroedinger equation for bound states with Mathematica 3.0 (arXiv link). – Jens Jan 31 '16 at 21:21
• Your code doesn't run because Alpha is undefined. Where did you get the values in the last table? – Jens Jan 31 '16 at 21:34
• I was initially confused about your use of Neumann boundary conditions. I think you must mean something else, because by comparing with your desired eigenvalues I concluded that you're really solving the radial equation for the reduced wave function, where no Neumann boundary condition is needed. The biggest problem is how to deal with the conditions of vanishing amplitude at infinity. – Jens Feb 1 '16 at 1:26
• SORRY,my bad,alpha is a pre-defined value,actually it's the schroedinger constant.Here it's value is 2..I'm terribly sorry that you have to find out this on your own.And thank you very much for your tutorial in NDEigensystem. – Turgon Feb 1 '16 at 12:49
• @Jens OMG, there are Goto[]s in that Arxiv paper's code!!! For a moment I thought I was reading a BASIC program. – Peltio Feb 8 '16 at 17:11

## 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;

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} =
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;
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} =
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.

• Thank you so much.So my problem is with the shift.In your code,without a shift,the eigenvalues are just like what I solved before,with little difference.BUT why do I need a shift?You said "sort" the eigenvalues,I don't know if I'm getting this right,is it because there are lots of redundant eigenvalues and the program needs to sort them and put the biggest as well as the most accurate ones ahead?And what if we use methods other than "Arnoldi"? – Turgon Feb 1 '16 at 14:18
• And about the accuracy,it seems decreased with higher energy level?The first few are nearly perfect for me,then it just goes very wrong.What do you think if we add f'[d]==-10^-6` as a boundary condition?The function would converge faster,I think.And this is what I need a Neumann condition for. – Turgon Feb 1 '16 at 14:45
• @user21 Standing on the shoulders of giants... – Jens Feb 1 '16 at 20:40
• I think I am going to print that one out - for bad days... Thanks. No, seriously if you have something I should look at concerning NDEigensystem, please let me know. Ideas for improvement, comments etc... – user21 Feb 1 '16 at 21:12
• The higher the excitation, the more nodes the wavefunction has. These are now squeezed into a fixed interval, so you need increased spatial resolution to reduce the error for high-lying states. Ignore the "brain" remark, it was a joke about something I don't use very often (but biologists claim everybody has one - it's just a theory, though). – Jens Feb 3 '16 at 17:17