Python runs in your browser. Heavy or infinite loops may freeze your browser tab. Use "Stop" if needed; for heavy jobs, run locally.

Program 1

Harmonic Oscillator

Given a one-dimensional quantum harmonic oscillator described by the potential:

\[ V(x) = \frac{1}{2} k x^2 \]

Solve the time-independent Schrödinger equation:

\[ -\frac{\hbar^2}{2m} \frac{d^2 \psi(x)}{dx^2} + V(x) \psi(x) = E \psi(x) \]

with the boundary conditions:

\[ \psi(-L) = 0 \quad \text{and} \quad \psi(L) = 0 \]

Solve 1-dimensional time-independent Schroedinger Equation (TISE) for Quantum Harmonic Oscillator (QHO) using shooting and odeint method to find energy eigenvalues and corresponding wavefunctions

import numpy as np
from scipy import integrate
from scipy.integrate import odeint, trapezoid
from scipy.optimize import bisect
import matplotlib.pyplot as plt

L, m, hbar, w, k, E0, En = 6, 1, 1, 1, 1, 0, 9
Erange = np.arange(E0, En, 0.1)
N = 1000
x, h = np.linspace(-L, L, N, retstep = True)
def V(x): return 0.5*k*x**2 # define potential   
psi = np.zeros(N)

u = [0, 0.1]  # Initial conditions for wavefunction and its derivative

def f(u, x, E):
    y, z = u
    f1 = z
    f2 = ((2 * m) / (hbar**2)) * (V(x) - E) * y
    return (f1, f2)

def psiboundval(E):
    sol = odeint(f, u, x, args=(E,))
    return sol[:, 0][-1]  # Boundary condition at r=L

def shoot(Erange):
    Y = np.array([psiboundval(E) for E in Erange])
    eigval = np.array([bisect(psiboundval, Erange[i], Erange[i + 1])
                       for i in np.where(np.diff(np.signbit(Y)))[0]])
    return eigval

# Find the energy eigenvalues
En = shoot(Erange)
print("Energy Eigen Values:", En)

# Plot the normalized wavefunctions for each energy eigenvalue
plt.figure(figsize=(10,15))
st = 0
# obtain the wave functions
for eigen in En:
    sol = odeint(f, u, x, args=(eigen,))
    psi = sol[:, 0]
    normpsi = psi/np.sqrt(trapezoid(psi*psi, x)) #Normalize the wave function
    plt.subplot(int(len(En) / 2) + 1, 2, st+1)
    plt.ylabel(r"$\Psi_"+str(st)+"(x)$", fontsize=16)
    plt.xlabel(r'$x$', fontsize=16)
    plt.axhline(y = eigen, color="black")
    plt.axvline(color="black")
    plt.xlim(-L,L)
    plt.ylim(-1, 10)
    plt.plot(x, eigen + normpsi, '-', color="red", label='$\Psi_'+str(st)+'(x)$') #plot wave function
    plt.plot(x, V(x))
    plt.legend()
    st = st + 1
plt.tight_layout()
plt.savefig('psi.png')
plt.show()

Output 1

Program 2

Harmonic Oscillator

Solve 1-dimensional time-independent Schroedinger Equation (TISE) for Quantum Harmonic Oscillator (QHO) using shooting and Numerov method

import numpy as np
from scipy import integrate
from scipy.integrate import odeint, trapezoid
from scipy.optimize import bisect
import matplotlib.pyplot as plt

L, m, hbar, w, k, E0, En = 6, 1, 1, 1, 1, 0, 9
E_range = np.arange(E0, En, 0.1)
N = 1000
x, h = np.linspace(-L, L, N, retstep = True)
V = 0.5*k*x**2 # define potential   
psi = np.zeros(N)

#  Integrate function with inital value psi0 & potential V. 
def numerov(E): 
    k2 = (2*m/hbar**2)*(E - V)  
    psi[0], psi[1] = 0, 0.05
    for i in range(1, N-1):
        psi[i+1] = (2*(1 - (5/12)*h**2*k2[i])*psi[i] - (1 + (1/12)*h**2*k2[i-1])*psi[i-1])/(1 + (1/12)*h**2*k2[i+1])
    return psi

# Define function to get right bound value of psi
def shoot_optim(E):
    psi = numerov(E)
    return psi[-1]

#define shooting function to get the energy eigen values
def shoot(Erange): 
    Y = np.array([shoot_optim(E) for E in Erange])
    eigval = np.array([bisect(shoot_optim, Erange[i], Erange[i + 1])
                       for i in np.where(np.diff(np.signbit(Y)))[0]])
    return eigval

En = shoot(E_range)
print("Energy Eigen Values:", En)
plt.figure(figsize=(10,15))
st = 0
# obtain the wave functions
for eigen in En:
    psi = numerov(eigen)
    normpsi = psi/np.sqrt(trapezoid(psi*psi, x)) #Normalize the wave function
    plt.subplot(int(len(En) / 2) + 1, 2, st+1)
    plt.ylabel(r"$\Psi_"+str(st)+"(x)$", fontsize=16)
    plt.xlabel(r'$x$', fontsize=16)
    plt.title('Energy level: n= '+str(st)+', Energy eigen value ='+str(eigen) )
    plt.axhline(y = eigen, color="black")
    plt.axvline(color="black")
    plt.xlim(-L,L)
    plt.ylim(-1, 10)
    plt.plot(x, eigen + normpsi, '-', color="red", label='$\Psi_'+str(st)+'(x)$') #plot wave function
    plt.plot(x, V)
    plt.legend()
    st = st + 1
plt.tight_layout()
plt.savefig('psi.png')
plt.show()

Output 2

Interactive Example Link

Click for more Interactive Examples

TISE Solution

Harmonic Oscillator

Harmonic Oscillator

Interactive Example Link