0
$\begingroup$

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[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}
  1. 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
Improve this question
$\endgroup$
15
  • $\begingroup$ What is the value of Hop0? $\endgroup$
    – Ulrich Neumann
    Dec 9, 2022 at 10:36
  • $\begingroup$ @Ulrich Neumann, thanks, I fixed it, it was a typo $\endgroup$
    – Mam Mam
    Dec 9, 2022 at 10:40
  • $\begingroup$ threes error in NIntegrate try to resolve this error firstly ! Think of alternative ways. your question posted 3 days ago. it still NIntegrate issue $\endgroup$
    – nufaie
    Dec 9, 2022 at 13:49
  • $\begingroup$ @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? $\endgroup$
    – Mam Mam
    Dec 9, 2022 at 14:19
  • $\begingroup$ 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. $\endgroup$
    – Roman
    Dec 9, 2022 at 14:19

1 Answer 1

Reset to default
1
$\begingroup$

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}    *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – Mam Mam
    Dec 9, 2022 at 16:01
  • $\begingroup$ 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. $\endgroup$
    – Roman
    Dec 9, 2022 at 17:33
  • $\begingroup$ 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? $\endgroup$
    – Mam Mam
    Dec 9, 2022 at 18:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.