# ground state from the Schroedinger equation with a central potential what happens to the origin

I have code that attempts to implement a solution to the Schrödinger equation where there is a central potential (more or less im thinking of hydrogen), in 1-D using the numerov method to construct the wave function. It does not take care of the states (N=1,..) as my goal is to understand the link between the theory and the software. Here are some results for some energies, whether or not they are eigenvalues, im not sure:

[6.054999999999915, 6.214999999999912, 6.534999999999905]


and a graph . The algorithm im trying to implement is noted here in section 6 PDF, but my potential is different. My problem is, the central potential has a 1/r at the origin, how would I deal with that? (I did a sidestep in the code by adding an arbitrary 0.00001 to eliminate 0).

There is a more general code in another question seen here, but I do not follow the shift from equation to 3d matrix.

Here is my code:

import numpy as np
import scipy
import matplotlib.pyplot as plt

h=np.float64(1)
m=60
c=m/2
estep = 1.0/(100)
delta = 0.0002
x=np.linspace(-1000,1000,num=m)
x=(x - x[m/2]+0.00001 )
v =  1.0/(x**2)
sig=None
E = []

def get_left(wf,k):
for i in range(2,int(c)):
wf[i] = 2*(( 1- 1/12* h**2* k[i-2]**2)*wf[i-2] - (1 + 1/12 *h**2* k[i-1]**2)*wf[i-1])/(1 + 1/12* h**2* k[i]**2)
wf[0:c] = wf[0:c]/(np.sqrt(np.sum(wf[0:c]**2)))
return wf

def get_right(wf,k):
mi = m-1
for i in range(2,int(c)):
wf[mi-i] = 2*(( 1- 1/12* h**2* k[mi-i+2]**2)*wf[mi-i+2] - (1 + 1/12 *h**2* k[mi-i+1]**2)*wf[mi-i+1])/(1 + 1/12* h**2*k[mi-i]**2)
wf[c:mi] = wf[c:mi]/(np.sqrt(np.sum(wf[c:mi]**2)))
return wf

def root_search(epsilonu, epsilonl):
en = None
it =0
while it < 10000:
phil =  np.zeros(m)
phir =  np.zeros(m)
phil =0
phil = 0.001
phir[m-1]=0
phir[m-2]=0.001
k = (-v + (epsilonu+epsilonl)/2)
phil = get_left(phil,k)
phir = get_right(phir,k)
diffl = np.diff(phil)
diffr = np.diff(phir)
erri = diffl[m/2]- diffr[m/2]
if erri < delta:
plt.plot(x, (phil+phir)**2 )
plt.savefig( str('sho'+'.png' ))
return (epsilonu+epsilonl)/2
if erri < 0 :
epsilonl = (epsilonu+epsilonl)/2
if erri > 0:
epsilonu =epsilonl
epsilonl = (epsilonu+epsilonl)/2
it+=1
return None

def start():
epsilonl=0.01
eps = epsilonl
epsilonu=None
iteration=0
sig = None
while iteration < 100000:
phil =  np.zeros(m)
phir =  np.zeros(m)
phil =0
phil = 0.001
phir[m-1]=0
phir[m-2]=0.001
k = (-v + eps)
phil = get_left(phil,k)
phir = get_right(phir,k)
diffl = np.diff(phil)
diffr = np.diff(phir)
err = diffl[m/2]- diffr[m/2]
ssig = 1 if err>0 else -1
if sig:
if sig != ssig:
epsilonu= eps
E.append(root_search(epsilonu, epsilonl))
print(E)
sig = ssig
epsilonl = eps
eps = eps+estep
iteration+=1

if __name__ == "__main__":
start()


The theory: Im attempting to solve the 1-D Schroedinger equation, $$\left( -\frac{h^2}{2m}\frac{d^2}{dx^2} + V\left(x\right)\right) \phi\left(x\right) = \epsilon \phi\left( x\right)$$ with V as a central potential (coulombic force, inverse square relation), the final equation looks like (with proper choice of units to give q^2=1, h-1..., energy in hartree) $$\frac{d^2}{dx^2}\phi\left(x\right) + 2 \left(\epsilon -\frac{1}{x^2} \right) \phi\left(x\right) =0$$

• If you're thinking about the hydrogen atom, then there shouldn't be a singularity at the origin, because you should do a separation of variables on the 3d PDE, instead of plugging in $V=1/r$ into the 1d PDE, and then there is no singularity at $r=0$ (also, $r\geq0$ then). Otherwise, how to handle a singular potential is more of a mathematics and physics question, rather than numerics: the wavefunction has to vanish so that $V\psi \neq \infty$ is consistent with the equation, and the probability current has to be conserved — that's really a mathematics question. – Kirill Dec 31 '15 at 8:44
• See, for example physics.stackexchange.com/q/69889/3696 – Kirill Dec 31 '15 at 8:48
• I don't think it's only a mathematics question -- properly handling functions with singularities is also a computational challenge. If possible, you want to rewrite the equation to remove the singularity (as you suggested) or otherwise take an analytic limit, I guess. – Jannis Teunissen Dec 31 '15 at 15:27
• I don't know what you mean by: but I do not follow the shift from equation to 3d matrix. If you mean how to translate the differential equation into a system of equations, it was done using Finite Differences. – nicoguaro Jan 2 '16 at 21:08
• If it is really a 1d problem, as far as I remember, already the 1/x potential leads to a diverging eigenstate, and then even more so the 1/x^2. As a workaround people sometimes introduce a threshold parameter, e.g. 1/(x^2-a)^0.5. But, if your equation is coming from a 3d background, consider Nicoguaro's points. – davidhigh Jun 18 '19 at 19:55

As mentioned by @Kirill, the differential equation for the hydrogen is slightly different. After the separation of variables, you end up with

$$\left[ - \frac{\hbar^2}{2\mu} \left({1 \over r^2}{\partial \over \partial r}\left(r^2 {\partial R(r)\over \partial r}\right) - {l(l+1)R(r)\over r^2} \right) + V(r)R(r) \right]= E R(r),$$

You can rewrite the equation as

$$\left[ -\frac{1}{2} \frac{d^2}{dy^2} + \frac{1}{2} \frac{l(l+1)}{y^2} - \frac{1}{y}\right] u_l = W u_l$$

The next code (based on my previous post) computes the eigenvalues/eigenvectors for this equation. See that I added a small value to the denominators, so they don't vanish at zero.

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

L = 50*np.pi
N = 1001
xa = 0
xb = L
x = np.linspace(xa, xb, N)
dx = x - x

l = 1
k = 0
T = -0.5*diags([-2., 1., 1.], [0, -1, 1], shape=(N, N))/dx**2
U_vec = 0.5*l*(l + 1)/(x**2 + 1e-6) - 1/(np.abs(x) + 1e-6)
U = diags([U_vec], )

H = T + U

vals, vecs = eigsh(H, which='SA')

print(np.round(vals, 6))
print(np.round([-1/(2*n**2) for n in range(k + l + 1, k + l + 7)], 6))

for k in range(5):
vec = vecs[:, k]
mag = np.sqrt(np.dot(vecs[:, k],vecs[:, k]))
vec = vec/mag
plt.plot(x, vec, label=r"$$n=%i$$"% (k+1))

plt.xlabel(r"$$x$$")
plt.ylabel(r"$$\psi(x)$$")
plt.xlim(xa, xb)
plt.legend()
plt.savefig("eigenvectors.png", dpi=600)
plt.show()


That gives as solution:

[-0.125106 -0.0556   -0.031272 -0.020012 -0.013896 -0.010209]
[-0.125    -0.055556 -0.03125  -0.02     -0.013889 -0.010204] For your repulsive potential (different sign than the Laplacian), I used the same code

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

L = 50*np.pi
N = 1001
xa = 0
xb = L
x = np.linspace(xa, xb, N)
dx = x - x

T = -0.5*diags([-2., 1., 1.], [0, -1, 1], shape=(N, N))/dx**2
U_vec = 1/(x**2 + 1e-6)
U = diags([U_vec], )

H = T + U

vals, vecs = eigsh(H, which='SA')

print(np.round(vals, 6))

for k in range(5):
vec = vecs[:, k]
mag = np.sqrt(np.dot(vecs[:, k],vecs[:, k]))
vec = vec/mag
plt.plot(x, vec, label=r"$$n=%i$$"% (k+1))

plt.xlabel(r"$$x$$")
plt.ylabel(r"$$\psi(x)$$")
plt.xlim(xa, xb)
plt.legend()
plt.savefig("eigenvecs.png", dpi=600)
plt.show()


[ 0.000408  0.001207  0.002405  0.004001  0.005997  0.008392] 