Numerical integration of Schroedinger equation

picb_integrate is specifically for particle-in-a-complicated-box, i.e., a one-sided finite square well

$$ -\frac{\hbar^2}{2m}\frac{d^2 \psi}{dx^2} + U(x) = E\psi $$

Using dimensionless length

• $x^\prime \equiv x/\sqrt{\hbar/2m\omega}$, and dimensionless energies
• $E^\prime = 2E/\hbar\omega$, Schroedinger equation becomes

$$ -\frac{2m\omega}{\hbar}\frac{d^2 \psi}{{dx^\prime}^2} + U\psi = E\psi, $$$$ \frac{d^2 \psi}{{dx^\prime}^2} = -\left(E^\prime - {x^\prime}^2\right) \psi. $$

The integration function uses a second order approximation for second derivative.

Developed for PHYS 212, March 2016 Modified for Jupyter notebook, Spring 2018
Slight updates December 2020

Marty Ligare

import numpy as np

import matplotlib as mpl
import matplotlib.pyplot as plt

# Following is an Ipython magic command that puts figures in the  notebook.
%matplotlib notebook

# M.L. modification of matplotlib defaults
# Changes can also be put in matplotlibrc file, 
# or effected using mpl.rcParams[]
mpl.style.use('classic')
plt.rc('figure', figsize = (6, 4.5))    # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size

def u(x):
    return x**2

def integrate(x):
    y = np.zeros(len(x))      # Initialize array for values of psi
    y[0] = 1
    y[1] = 1
    for i in range(1,n-1):   # Update values of y[2], y[3], y[4], etc.
        y[i+1] = 2*y[i] - y[i-1] - dx**2*(e-u(x[i]))*y[i]  
    return y  

"Guess" for e (energy)

e = 9.00001  # "Guess" for value of dimensionless energy

lep = 0                       # Left end-point 
rep = lep + 6                 # Right end-point 
n = 1000001                   # Number of points (number of intervals = np - 1)
dx = (rep-lep)*1./(n-1)      # Distance between points

x = np.linspace(lep, rep, n) # Create array of x-values
y = integrate(x)              # Create values for wavefunction

plt.figure()

plt.xlabel('$x^\prime$')      # Label for horizontal axis
plt.ylabel("$\psi$")    # Label for vertical axis
plt.title("simple harmonic oscillator")
plt.grid(True)
plt.axhline(0,color='green')  # Makes solid green x-axis
plt.ylim(-2,2)
plt.plot(x,y);

The "guess" of $E^\prime = 9.00001$ gives decent as $x^\prime$ gets large. Recall

$$ E = \frac{1}{2}\hbar\omega E^\prime = \frac{1}{2}\hbar\omega \times 9 = 4.5\hbar\omega $$

which is the corrrect energy for the illustrated $n=4$ state.

Version information

• %version_information is an IPython magic extension for showing version information for dependency modules in a notebook;

• See https://github.com/jrjohansson/version_information

• %version_information is available on Bucknell computers on the linux network. You can easily install it on any computer.

%load_ext version_information

version_information numpy, matplotlib

Software	Version
Python	3.7.7 64bit [GCC 7.3.0]
IPython	7.16.1
OS	Linux 4.9.0 9 amd64 x86_64 with debian 9.13
numpy	1.18.5
matplotlib	3.2.2
Mon Dec 21 11:20:37 2020 EST