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

TISE by Finite Difference Method

Solve the time-independent Schrödinger equation in one dimension:

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

where: $ \psi(x) $ is the wavefunction, $ V(x) $ is the potential energy, $ E $ is the energy eigenvalue.

Consider potentials:

Harmonic oscillator: $ V(x) = \frac{1}{2} m \omega^2 x^2. $

The objective is to compute the energy eigenvalues $ E_n $ and corresponding normalized wavefunctions $ \psi_n(x) $, and to calculate expectation values $ \langle x \rangle $ and $ \langle x^2 \rangle $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid
from numpy.linalg import eigh

m, hbar = 1, 1
def LHO(N, L, omega=1):
    x = np.linspace(-L/2, L/2, N)
    V = np.zeros(N)
    V = 0.5 * m * omega**2 * x**2
    return x, V

L, N = 6, 1000
x, V = LHO(N, L)   # Potential
dx = L / (N - 1)
# Kinetic energy term (second derivative)
T = (np.eye(N)*(-2) + np.eye(N, k=1) + np.eye(N, k=-1)) / dx**2
H = -(hbar**2 / (2*m)) * T + np.diag(V)
eigenvalues, eigenvectors = eigh(H)
    
# Normalize each eigenfunction
for i in range(len(eigenvectors)):
    norm = np.sqrt(trapezoid(np.abs(eigenvectors[:, i])**2, x))
    eigenvectors[:, i] /= norm

# Print first 3 eigenvalues
print(f"First 3 eigenvalues: {eigenvalues[:3]}")

# Plot the first few wavefunctions
plt.plot(x, V / np.abs(V).max(), 'k--', label="Scaled Potential")
plt.plot(x, np.abs(eigenvectors[:, 0])**2, label=r'$\psi_0(x)^2$', color='blue')
plt.plot(x, np.abs(eigenvectors[:, 1])**2, label=r'$\psi_1(x)^2$', color='red')
plt.plot(x, np.abs(eigenvectors[:, 2])**2, label=r'$\psi_2(x)^2$', color='green')
plt.title("Wavefunctions and Potential (LHO)")
plt.xlabel("x")
plt.ylabel("Probability Density (normalized)")
plt.legend()
plt.grid(True)
plt.savefig("plot.png")
plt.show()

x_expect = trapezoid(eigenvectors[:, 0].conj()*x*eigenvectors[:, 0], x)
x2_expect = trapezoid(eigenvectors[:, 0].conj()*x**2*eigenvectors[:, 0], x)

print(f"Expectation value of x for ground state = {x_expect:.4f}")
print(f"Expectation value of x2 for ground state = {x2_expect:.4f}")

Output 1

Finite Difference Approximation

The TISE is:

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

The second derivative is approximated using central differences:

\[ \frac{d^2\psi}{dx^2}\Big|_{x_i} \approx \frac{\psi_{i+1} - 2\psi_i + \psi_{i-1}}{(\Delta x)^2} \] Thus, we get \[\small -\frac{\hbar^2}{2m(\Delta x)^2}[\psi_{i+1}+(-2+V_i)\psi_i+\psi_{i-1}] = E\psi_i \]

where $ \Delta x = L/(N - 1) $ is the grid spacing.

Hamiltonian in Matrix Form

Discretizing the spatial domain into $ N $ points, the wavefunction becomes a column vector:

\[ \vec{\psi} = \begin{bmatrix} \psi_1 \\[4pt] \psi_2 \\[2pt] \vdots \\[2pt] \psi_N \end{bmatrix} \]

The kinetic energy operator using finite differences is represented by the tridiagonal matrix $ T $:

\[\small T = \frac{1}{(\Delta x)^2} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \cdots & 0 \\ 0 & 1 & -2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 1 & -2 \end{bmatrix} \]

The potential energy operator becomes a diagonal matrix:

\[\small V = \begin{bmatrix} V_1 & 0 & \cdots & 0 \\ 0 & V_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & V_N \end{bmatrix} \]

The Hamiltonian matrix is then:

\[ H = -\frac{\hbar^2}{2m} T + V \]

The time-independent Schrödinger equation in matrix form is:

\[ H \vec{\psi} = E \vec{\psi} \]

or explicitly, component-wise:

\[\tiny \begin{bmatrix} H_{11} & H_{12} & 0 & \cdots & 0 \\ H_{21} & H_{22} & H_{23} & \cdots & 0 \\ 0 & H_{32} & H_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & H_{N-1,N} \\ 0 & 0 & 0 & H_{N,N-1} & H_{NN} \end{bmatrix} \begin{bmatrix} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_N \end{bmatrix} = E \begin{bmatrix} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_N \end{bmatrix} \]

where

\[\small H_{ii} = \frac{2\hbar^2}{2m(\Delta x)^2} + V_i, \quad H_{i,i\pm1} = -\frac{\hbar^2}{2m(\Delta x)^2}, \]

and all other elements are zero.

Solving this eigenvalue problem yields the allowed energies $ E $ and the corresponding normalized wavefunctions $ \vec{\psi} $.

Expectation Values

Expectation values are computed numerically as:

\[ \langle x \rangle = \int x |\psi(x)|^2 \, dx \] \[ \langle x^2 \rangle = \int x^2 |\psi(x)|^2 \, dx \]

Algorithm

Define the potential $ V(x) $ on a discrete grid of $ N $ points over domain $ L $.
Construct the kinetic energy matrix $ T $ using finite differences (tridiagonal matrix).
Construct the Hamiltonian matrix: $ H = -\frac{\hbar^2}{2m}T + \text{diag}(V) $.
Solve the eigenvalue problem $ H\vec{\psi} = E\vec{\psi} $ using numpy.linalg.eigh.
Normalize the wavefunctions.
Plot probability densities $ |\psi_n(x)|^2 $.
Compute expectation values $ \langle x \rangle $ and $ \langle x^2 \rangle $ using the trapezoidal rule.

Program 2

TISE by Finite Difference Method

Solve the time-independent Schrödinger equation in one dimension for the harmonic oscillator potential $ V(x) = \frac{1}{2} m \omega^2 x^2 $:

by approximate the second derivative using the finite difference method:

\[ \frac{d^2\psi}{dx^2}\Big|_{x_i} \approx \frac{\psi_{i+1} - 2\psi_i + \psi_{i-1}}{(\Delta x)^2}, \]

Solve the eigenvalue equation:

\[ H\vec{\psi} = E\vec{\psi} \]

to obtain the energy eigenvalues $ E_n $ and normalized eigenfunctions $ \psi_n(x) $.

Normalize each eigenfunction numerically:

\[ \int |\psi_n(x)|^2 dx = 1 \]

Plot the potential $ V(x) $ and the first three probability densities.

print:

The first three eigenvalues $ E_0, E_1, E_2 $
Normalization checks $ \langle \psi_n | \psi_n \rangle $
Orthogonality checks $ \langle \psi_0 | \psi_1 \rangle $ and $ \langle \psi_0 | \psi_2 \rangle $

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid
from numpy.linalg import eigh

m, hbar, omega = 1, 1, 1
def LHO(N, L):
    x = np.linspace(-L/2, L/2, N)
    V = np.zeros(N)
    V = 0.5 * m * omega**2 * x**2
    return x, V

L, N = 12, 1000
x, V = LHO(N, L)   # Potential
dx = L / (N - 1)
# Kinetic energy term (second derivative)
T = (np.eye(N)*(-2) + np.eye(N, k=1) + np.eye(N, k=-1)) / dx**2
H = -(hbar**2 / (2*m)) * T + np.diag(V)
eigenvalues, eigenvectors = eigh(H)
    
# Normalize each eigenfunction
for i in range(len(eigenvectors)):
    norm = np.sqrt(trapezoid(np.abs(eigenvectors[:, i])**2, x))
    eigenvectors[:, i] /= norm

# Print first 3 eigenvalues
print(f"First 3 eigenvalues: {eigenvalues[:3]}")

# Plot the first three wavefunctions
plt.plot(x, V , 'k--', label="Potential")
plt.plot(x, eigenvalues[0] + np.abs(eigenvectors[:, 0])**2, label=r'$\psi_0(x)^2$', color='blue')
plt.plot(x, eigenvalues[1] + np.abs(eigenvectors[:, 1])**2, label=r'$\psi_1(x)^2$', color='red')
plt.plot(x, eigenvalues[2] + np.abs(eigenvectors[:, 2])**2, label=r'$\psi_2(x)^2$', color='green')
plt.title("Wavefunctions and Potential (LHO)")
plt.xlabel("x")
plt.ylabel("Scaled Probability Density (normalized)")
plt.xlim(-L/4, L/4)
plt.ylim(-0.2, np.max(eigenvalues[:3]) + 0.5)
plt.legend()
plt.grid(True)
plt.savefig("plot.png")
plt.show()

#Check Orthonormal conditions
norm00 = trapezoid(np.abs(eigenvectors[:, 0])**2, x)
norm11 = trapezoid(np.abs(eigenvectors[:, 1])**2, x)
norm01 = trapezoid(eigenvectors[:,0].conj()*eigenvectors[:, 1], x)
norm02 = trapezoid(eigenvectors[:,0].conj()*eigenvectors[:, 2], x)

print(f"Normalization check for ground state = {norm00:.4f}")
print(f"Normalization check 1st excited state = {norm11:.4f}")
print(f"Orthogonal condition check between ground and 1st excited state = {norm01:.4f}")
print(f"Orthogonal condition check between ground and 2nd excited state = {norm02:.4f}")

Output 2

Program 3

TISE by Finite Difference Method

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import eigh
from scipy.integrate import trapezoid
m, hbar = 1, 1
def Hatom(N, L, e=1):
    x = np.linspace(0.01, L, N)
    V = np.zeros(N)
    V = -e**2 / x
    return x, V

L, N = 20, 1000
x, V = Hatom(N, L)   # Potential
dx = L / (N - 1)
# Kinetic energy term (second derivative)
T = (np.eye(N)*(-2) + np.eye(N, k=1) + np.eye(N, k=-1)) / dx**2
H = -(hbar**2 / (2*m)) * T + np.diag(V)
eigenvalues, eigenvectors = eigh(H)
    
# Normalize each eigenfunction
for i in range(len(eigenvectors)):
    norm = np.sqrt(trapezoid(np.abs(eigenvectors[:, i])**2, x))
    eigenvectors[:, i] /= norm

# Print first 3 eigenvalues
print(f"First 3 eigenvalues: {eigenvalues[:3]}")

# Plot the first few wavefunctions
plt.plot(x, np.abs(eigenvectors[:, 0])**2, label=r'$\psi_0(x)^2$', color='blue')
plt.plot(x, np.abs(eigenvectors[:, 1])**2, label=r'$\psi_1(x)^2$', color='red')
plt.plot(x, np.abs(eigenvectors[:, 2])**2, label=r'$\psi_2(x)^2$', color='green')
plt.title("Wavefunctions and Potential (Hydrogen-like)")
plt.xlabel("x")
plt.ylabel("Probability Density (normalized)")
plt.legend()
plt.grid(True)
plt.savefig("plot.png")
plt.show()

# Expectation values
x_expect = trapezoid(eigenvectors[:,0].conj()*x*eigenvectors[:, 0], x)
x2_expect = trapezoid(eigenvectors[:,0].conj()*x**2*eigenvectors[:, 0], x)

print(f"< x > for ground state = {x_expect:.4f}")
print(f"< x2 > for ground state = {x2_expect:.4f}")

Output 3

TISE Solution

TISE by Finite Difference Method

Finite Difference Approximation

Hamiltonian in Matrix Form

Expectation Values

Algorithm

TISE by Finite Difference Method

TISE by Finite Difference Method