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?
- numerically
In[82]:= 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[84]:= (Take[Sort[vals], 5] - 0.5)
Out[84]= {-2.00001, -0.500001, -0.222223, -0.125, -0.0800001}
- matrix
In[26]:= 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[32]= {-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[33]:=
Emin = Min[Eeig]
Out[33]= -1.5584632981876349185146804146007
Hop0
? $\endgroup$NDEigensystem
approach you seem to be calculating the eigensystem ofH + Psi[r]*0.5
instead ofH
, 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. $\endgroup$