1
$\begingroup$

this is my first question on here.

I'm trying to solve the Schrödinger equation for the hydrogen atom in the following form numerically:

$$\left[-\frac{\hbar^2}{2m}\frac{d^2}{dr^2}+V(r)+\frac{\hbar^2l(l+1)}{2mr^2}\right]R(r)=ER(r).$$

With the Coulomb-Potential $$V(r) = -\frac{e^2}{4\pi\epsilon_0r^2}.$$ By introducing the dimensionless variables $x = r/a_0$ and $\epsilon = E/E_0$, $E_0=1\text{Ry}$, the Schrödinger equation can be written as follows:

$$\left[-\frac{d^2}{dx^2}-\frac{2}{x}+\frac{l(l+1)}{x^2}\right]R(x)=\epsilon R(x)$$

I used the method of Leandro M. described in this post to solve the Schrödinger equation by discretizing the derivative of $R$ and transforming it into a Matrix equation.

I implemented this solution in Python and I think it worked, the only problem I have is that I don't know how to get from the Eigenvalues and -vectors of the resulting Matrix to the Energies and the corresponding Wavefunctions of the Hydrogen Atom for $n=1,2,3,...$ and $l = 0$.

Here is my code:

import numpy as np
from numpy.linalg import eig
from matplotlib import pyplot as plt

d = 0.05
# set values for r
rsteps = 3000
rmax = 100
r = np.linspace(1e-7, rmax, rsteps)

# create first matrix of schrodinger eq corresponding to the derivative with the following shape:
# (-2  1                   )
# ( 1 -2  1                )
# (    1 -2  1             )  * (1/(d**2))
# (         ...            )
# (                1  -2  1)
# (                    1 -2)


m = np.empty((len(r), len(r)))

x = -1 / d ** 2

m[0, 0] = -2 * x
m[0, 1] = 1 * x
m[len(r) - 1, len(r) - 1] = -2 * x
m[len(r) - 1, len(r) - 2] = 1 * x

for i in range(1, len(r) - 1):
    m[i, i - 1] = 1 * x
    m[i, i] = -2 * x
    m[i, i + 1] = 1 * x

# create matrix corresponding to the potential with the following shape:
# (V1                 )
# (   V2              )
# (      V3           )
# (       ...         )
# (            VN-1   )
# (                 VN)

vdiag = []
l = 0

for i in range(0, len(r) - 1):
    vdiag.append(-2 / r[i] + l * (l + 1) / (r[i] ** 2))

v = np.empty((len(r), len(r)))
np.fill_diagonal(v, vdiag)

# add matrices to get eigenvalue equation: H*R(x) = E*R(x)
H = np.empty((len(r), len(r)))

for i in range(0, len(v[0] - 1)):
    for j in range(0, len(v[1] - 1)):
        H[i, j] = m[i, j] + v[i, j]

# setting boundary conditions R_1 = R_N = 0
H[:, 0] = H[:, len(H[1]) - 1] = 0

# determine eigenvalues and eigenvectors
energies, wavefcts = eig(H)
$\endgroup$
6
  • 5
    $\begingroup$ I’m voting to close this question because it is about debugging code. $\endgroup$
    – Jon Custer
    Nov 18 '21 at 16:35
  • 1
    $\begingroup$ I don't want you to debug my code, I just want to know the last step to get from where my code ends to the actual Eigenergies and Wavefunctions. $\endgroup$ Nov 18 '21 at 16:37
  • 2
    $\begingroup$ What do you mean "actual"? Are they not stored in 'energies' and 'wavefcts' after the last line...? $\endgroup$ Nov 18 '21 at 16:40
  • 1
    $\begingroup$ you just need to reconvert the eigenvalues that you governed in the last line into the unit of energy by doing the transformation of the SE backwards. $\endgroup$
    – franz
    Nov 18 '21 at 16:56
  • 1
    $\begingroup$ Even if I'd do that, the resulting Energies would be way too big, so I think there is still a step missing to get the correct Eigenenergies. $\endgroup$ Nov 18 '21 at 17:34