CODEWITHSUNDEEP
pythonnumpy
Evan Walker
Evan Walker

Reputation: 25

Solving Shrodinger's equation for a particle in a harmonic potential well

Hello (this is my first time posting in stack overflow), I am trying the calculate the first 3 energy levels of a particle in a harmonic potential using the shooter method

The code is adapted from a script in Computational Physics by Mark Newman, this script calculated the ground state for a particle in a box.

here is the link http://www-personal.umich.edu/~mejn/cp/programs/squarewell.py

here is my adapted code

import numpy as np

from scipy.constants import m_e,hbar,elementary_charge

#define constants

vo = 50
a = 1e-11

#define limits of integration for adaptive runge-kutta method

xi = -10*a
xf = 10*a
N = 1000
dx = (xf-xi)/N

def V(x):#define harmonic potential well
    return vo*((x**2)/(a**2))

def f(r,x,E): #schrodinger's equation
    psi = r[0]
    phi = r[1]
    fpsi = psi
    fphi = (2*m_e/(hbar**2))*(V(x)-E)*psi
    return np.array([fpsi,fphi],float)

def solve(E): #calculates wave function for an energy E
    psi = 0.0
    phi = 1.0
    r = np.array([psi,phi],float)
    for x in np.arange(xi,xf,dx): #adaptive runge-kutta method
        k1 = dx*f(r,x,E)
        k2 = dx*f(r+0.5*k1,x+0.5*dx,E)
        k3 = dx*f(r+0.5*k2,x+0.5*dx,E)
        k4 = dx*f(r+k3,x+dx,E)
        r += (k1+2*k2+2*k3+k4)/6
    return r[0]

#finds the energy using secant method

E1 = 0.0
E2 = elementary_charge
psi2 = solve(E1)
target = elementary_charge/1000

while abs(E1-E2)>target:
    psi1,psi2 = psi2,solve(E2)
    E1,E2 = E2,E2-psi2*(E2-E1)/(psi2-psi1)
    print (E2/elementary_charge)

when run I get this error

RuntimeWarning: invalid value encountered in double_scalars

E1,E2 = E2,E2-psi2*(E2-E1)/(psi2-psi1)

which I think means that psi2 and psi1 are too close together but I am not quite sure how to fix this

Upvotes: 1

Views: 775

Answers (1)

Rachayita Giri
Rachayita Giri

Reputation: 487

Yep! You are right about the values being too close to each other. Your code returns a nan. It is because of the division by zero.

I would suggest using a correction factor. Something like max(delta, (psi2-psi1)) in the denominator where delta can still be a very small value but it will prevent the division by zero.

Upvotes: 2

Related Questions