In our assignment, we are required to find the energies of the ground state and the first two excited states of the Schrödinger equation in a harmonic potential:
$$V = \frac{50 x^2}{(10^{-11})^2}\, .$$
I have the Python code for solving Schrödinger equation using secant Runge Kutta method:
import numpy as np
m = 9.1094e-31
hbar = 1.0546e-34
e = 1.6022e-19
L = 5.2918e-11
N = 1000
h = L/N
def V1(x):
return 0.0
#return (50*x**2)/(10**-11)**2
def V2(x):
return 0.0
#return (50*x**4)/(10**-11)**4
def f(r, x, E, V):
psi = r[0]
phi = r[1]
fpsi = phi
fphi = (2*m/hbar**2)*(V(x)-E)*psi
return np.array([fpsi,fphi], float)
def solve(E, V):
psi = 0.0
phi = 1.0
r = np.array([psi, phi], float)
for x in np.arange(0,L,h):
k1 = h*f(r, x, E, V)
k2 = h*f(r+0.5*k1, x+0.5*h, E, V)
k3 = h*f(r+0.5*k2, x+0.5*h, E, V)
k4 = h*f(r+k3, x+h, E, V)
r += (k1 + 2*k2 + 2*k3 + k4)/6
return r[0]
print("Start program")
E1 = 0.0
E2 = e
psi2 = solve(E1, V1)
target = e/1000
while abs(E1-E2)>target:
psi1, psi2 = psi2, solve(E2, V1)
E1, E2 = E2, E2-psi2*(E2-E1)/(psi2-psi1)
print("For V(x) = V0x^2/a^2:")
print("E=", E2/e, "eV")
However, there are two problems that arise that causes me not being able to find the solution to this assignment. And I hope that someone can help me with it.
First of all, when the function V1(x) returns 0.0, the program accurately outputs the energy of ground state for V1(x) = 0, which is E= 134.28637169369105 eV. However, when we change to V1(x) = V = (50x^2)/(10^-11)^2, which is what our assignment requires us to find, the program creates the following error:
C:/Users/wormw/Desktop/PHYS 3142 - Computational Methods in Physics/Assignment 5/Assignment 5.py:29: RuntimeWarning: overflow encountered in double_scalars
return np.array([fpsi,fphi], float)
C:/Users/wormw/Desktop/PHYS 3142 - Computational Methods in Physics/Assignment 5/Assignment 5.py:54: RuntimeWarning: invalid value encountered in double_scalars
And the output of E is E = nan eV. So, I have no idea why this would happen and how to fix it.
Secondly, the assignment also requires us to find the first two excited state, but I have no idea how to modify the code to obtain the eigenvalue for getting the excited states
Thirdly, the assignment also states that Try using x = -10a to +10a, with the wave function psi = 0 at both boundaries. Should I change the value of L at the start?
