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?

  • $\begingroup$ I suggest you change your system of units first. Try the one described in this post. $\endgroup$
    – nicoguaro
    Apr 16 '19 at 19:39
  • $\begingroup$ I agree, introduce a new variable x' = cx. Here, c can absorb the factor 50 and the denominator. Then do your numerics on x'. $\endgroup$
    – MPIchael
    Apr 17 '19 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.