# Why do different calculations give different eigenvalues?

I find the eigenvalues of the Hamiltonian in two ways - numerically and matrix . Why do different calculations give different eigenvalues?
The Hamiltonian is H = -1/2*Laplacian -(2/r)
It can be seen that the minimum energy calculated numerically is -2.00001, and matrix-wise is -1.55846. Why is this happening?

1. numerically
In:= H =
Simplify[-1/2*
Laplacian[Psi[r], {r, \[Theta], \[Phi]}, "Spherical"] -
Psi[r]*2/r];

{vals, funs} =
NDEigensystem[{H + Psi[r]*0.5}, Psi[r], {r, 0, 100}, 100,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.015}}},
"Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}];

In:= (Take[Sort[vals], 5] - 0.5)

Out= {-2.00001, -0.500001, -0.222223, -0.125, -0.0800001}


1. matrix
In:= ClearAll["Global*"]
Psi1[r_, n1_] = (2 E^(-(r/n1)) Sqrt[n1!/(-1 + n1)!] Hypergeometric1F1[
1 - n1, 2, (2 r)/n1])/n1^2;

Psi2[r_, n2_] = (2 E^(-(r/n2)) Sqrt[n2!/(-1 + n2)!] Hypergeometric1F1[
1 - n2, 2, (2 r)/n2])/n2^2;

(*kinetic energy*)

K[r_, n1_, n2_] =
FullSimplify[
Psi2[r, n2]*
Laplacian[Psi1[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];

(*potential energy*)

P[r_] = -2/r;

(*calculation of matrix elements,I am using 20 basis functions*)

EE = Table[-1/2*
NIntegrate[K[r, n1, n2]*r^2, {r, 0, \[Infinity]},
MaxRecursion -> 15, WorkingPrecision -> 32] +
NIntegrate[Psi2[r, n2]*P[r]*Psi1[r, n1]*r^2, {r, 0, \[Infinity]},
MaxRecursion -> 15, WorkingPrecision -> 32], {n1, 1, 20}, {n2,
1, 20}] // Chop;

(*eigenvalues*)

Eeig = Eigenvalues[EE]

Out= {-1.5584632981876349185146804146007, \
-0.3623635914395303877441420153245, \
-0.15669498305937198905543539839599, \
-0.08683775701950283493030149036754, \
-0.05501734483287334369748266587302, \
-0.03790244105929780721020220539079, \
-0.027651029147976237708049044437355, \
-0.021028273261568554013526089309235, \
-0.016502425356337140815385665090984, \
-0.013271998387649143824131249079151, \
-0.010884520722145291443991307150682, \
-0.009068777306934890312009174759412, \
-0.007654074742951805746421678011622, \
-0.006528495175126375499971622294854, \
-0.005615980338382413802381791780753, \
-0.004863061581404067179392515633070, \
-0.004230715950120777335031723062825, \
-0.003688843483376281574694297474038, \
-0.0032113000791884747824869858880628, \
-0.0027659547381622397716009741638790}

In:=
Emin = Min[Eeig]

Out= -1.5584632981876349185146804146007

• What is the value of Hop0? Dec 9, 2022 at 10:36
• @Ulrich Neumann, thanks, I fixed it, it was a typo Dec 9, 2022 at 10:40
• threes error in NIntegrate try to resolve this error firstly ! Think of alternative ways. your question posted 3 days ago. it still NIntegrate issue Dec 9, 2022 at 13:49
• @Alrubaie, I fixed the code as Daniel Huber advised me, now it is error-free. But this does not change the figures and the question remains the same, why are the minimum eigenvalues obtained by different methods for the same Hamiltonian different? Dec 9, 2022 at 14:19
• In the NDEigensystem approach you seem to be calculating the eigensystem of H + Psi[r]*0.5 instead of H, and then subtracting 0.5 from the eigenvalues. This is only equivalent if the normalization is equal to 1, which you don't enforce. To be exact, the Jacobian $r^2$ needs to be included when normalizing; but you are not telling the solver anything about this Jacobian. As a solution, you could try using $u(r)=r\cdot\psi(r)$ as radial wavefunctions, which can be normalized with a trivial Jacobian. Dec 9, 2022 at 14:19

A common trick in spherical coordinates is to define $$\psi(r)=\frac{u(r)}{r}$$ to simplify the normalization and Laplacian expressed in terms of $$u(r)$$:

$$\nabla^2\psi(r)=\psi''(r)+\frac{2}{r}\psi'(r) \qquad\leftrightarrow\qquad \nabla^2\frac{u(r)}{r}=\frac{u''(r)}{r}\\ \int_0^{\infty}\psi^2(r)r^2dr=1 \qquad\leftrightarrow\qquad \int_0^{\infty}u^2(r)dr=1$$

The Schrödinger equation then becomes $$-\frac12\nabla^2\psi(r)-\frac{2}{r}\psi(r)=E\psi(r) \qquad\leftrightarrow\qquad -\frac12u''(r)-\frac{2}{r}u(r)=E u(r)$$ which is much easier to deal with (looking like pseudo-straight coordinates with trivial Jacobian).

H = -1/2 u''[r] - 2/r u[r];

{vals, funs} = NDEigensystem[H, u[r], {r, 0, 100}, 100,
Method ->
{"SpatialDiscretization" -> {"FiniteElement",
{"MeshOptions" -> {"MaxCellMeasure" -> 0.015}}},
"Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}] //
Transpose // Sort // Transpose;

vals[[;; 5]]
(*    {-1.46467, -0.424888, -0.199048, -0.11502, -0.0748264}    *)
`
• Thanks a lot. The actual values have become closer to those obtained by the matrix method, but still they differ quite a lot. The minimum energy obtained is numerically is -1.46467, and matrix is -1.55846, but these values must be the same. Why do they turn out different? Dec 9, 2022 at 16:01
• The eigenvalues depend strongly on the spatial discretization. I suppose it's very difficult to discretize a grid well in order to properly treat a singular potential. Dec 9, 2022 at 17:33
• I'm new to this and some of the terms are a little tricky for me. By discretization you mean the area of changing r and MaxCellMeasure? Dec 9, 2022 at 18:04