TISE Solution - Harmonic Oscillator

Program 1

Harmonic Oscillator

Python 3.11 Ready

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()

Click "Run Code" to execute script and render figures.

Mathematical Problem Formulation

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

Program 2

Harmonic Oscillator

Python 3.11 Ready

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()

Click "Run Code" to execute script and render figures.

Mathematical Problem Formulation

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

Theoretical Background & Explanation

Interactive Example Link

Click for more Interactive Examples

Ch. 18.4 Harmonic Oscillator

Harmonic Oscillator

Mathematical Problem Formulation

Harmonic Oscillator

Mathematical Problem Formulation

Theoretical Background & Explanation

Interactive Example Link