TISE Solution: Harmonic Oscillator - Computational Physics Formulation & Simulation

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 & Computational Mechanics

Interactive Example Link

Click for more Interactive Examples

Theoretical Principles & Computational Physics Guide

Harmonic Oscillator: Differential Equations & Dynamical Simulation

Dr. Alorika Chatterjee & Dr. Anirban Shaw Science & Physics Learners Python 3 / NumPy / SciPy

1. Core Principles & Theoretical Framework

Dynamical Evolution: Solves differential equations governing kinematics, planetary orbits, electrical circuits, and wave propagation.

Step-by-Step Stepping: Advances continuous state vectors forward in discrete time intervals $\Delta t$ using slope estimates.

Conservation Laws: Evaluates preservation of physical invariants (total energy, angular momentum) across long integration horizons.

2. Mathematical Formulation & Algorithmic Pipeline

Vector initial value problem and numerical state advancement:

$$ \frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y}), \quad \mathbf{y}(t_0) = \mathbf{y}_0 \implies \mathbf{y}_{n+1} = \mathbf{y}_n + \Delta t \cdot \mathbf{\Phi}(t_n, \mathbf{y}_n, \Delta t) $$

3. Numerical Mechanics & Implementation Logic

Coupled Reduction: Convert $n$-th order ODEs into systems of $n$ coupled first-order equations.
Intermediate Slopes: Compute weighted intermediate derivatives ($k_1, k_2, k_3, k_4$ in RK4).
State Update: Advance state variables simultaneously and verify energy/phase invariants.

4. Relevant Computational Practice Exercises

Exercise 1

Solve the differential equation for Harmonic Oscillator over 1000 time steps and plot state coordinates vs time.

Exercise 2

Compare numerical solution stability for small vs large step sizes $h = \Delta t$.

Exercise 3

Plot phase space trajectory demonstrating physical stability or damping attractors.

Standard Literature & References:

John R. Taylor, "Classical Mechanics", University Science Books (2005).
Mark Newman, "Computational Physics", University of Michigan (2013).

Ch. 19.4 Harmonic Oscillator

Harmonic Oscillator

Mathematical Problem Formulation

Harmonic Oscillator

Mathematical Problem Formulation

Theoretical Background & Computational Mechanics

Interactive Example Link

1. Core Principles & Theoretical Framework

2. Mathematical Formulation & Algorithmic Pipeline

3. Numerical Mechanics & Implementation Logic

4. Relevant Computational Practice Exercises