# 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] =0
phil[1] = 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] =0
phil[1] = 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. Commented Dec 31, 2015 at 8:44
• See, for example physics.stackexchange.com/q/69889/3696 Commented Dec 31, 2015 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. Commented Dec 31, 2015 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. Commented Jan 2, 2016 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. Commented Jun 18, 2019 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[1] - x[0]

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], [0])

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[1] - x[0]

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], [0])

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]


Although, I am not sure what that means. Because it is not letting the electron to be in the origin pushing it away, and the domain should be really unbounded. It is similar to having a square potential with a round edge at the left. For the more realistic case, you can use nonuniform finite differences and increase them exponentially or absorbing boundary conditions.