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

Shooting Method to Solve TISE

Shooting Method to solve TISE

Consider a particle of mass $ m $ confined in a one-dimensional box of length $ L $. The potential energy $ V(x) $ inside the box is zero, while it is infinite outside the box. The goal is to solve the time-independent Schrödinger equation (TISE) for this system, which is given by:

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

where:

$ \hbar $ is the reduced Planck's constant,
$ E $ is the energy of the particle,
$ \Psi(x) $ is the wave function.

Implement a numerical solution using the shooting method to find the energy eigenvalues of the particle in the box.

'''Solve 1-dimensional time-independent Schroedinger Equation (TISE) for
particle in 1D box'''

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

def V(x):
    return 0

L, m, hbar, w, E0, En = 10, 1, 1, 1, 0.001, 5
Erange = np.arange(E0, En, 0.01)
x = np.linspace(0, L, 100)
u = [0, 0.01]
    
def f(u, x, E):
    y, z = u
    f1, f2 = z, ((2*m)/(hbar)**2)*(V(x) -E)*y
    return (f1, f2)

#Calculate psi(L)
def psiboundval(E):
    sol = odeint(f, u, x, args = (E, ))
    return sol[:, 0][-1]

#Identify eigenvalues (where psi(L) crosses zero)
Y = np.array([psiboundval(E) for E in Erange])
plt.plot(Erange, Y, label="$\Psi(L)$")
plt.axhline(y=0, color='red', linestyle='--', label="Boundary Condition $\Psi(L)$=0")
plt.ylabel(r"$\Psi(L)$")
plt.xlabel(r"$E$")
plt.title(r"$\Psi(L) vs. E$")
plt.legend()
plt.savefig("plot.png")
plt.show()

# Determine eigenvalues by Bisection method
En = []
for i in range(len(Y) - 1):
    if np.sign(Y[i]) != np.sign(Y[i + 1]):  # Check for sign change
        E_eigen = bisect(psiboundval, Erange[i], Erange[i + 1])
        En.append(E_eigen)
print("Energy Eigen Values:", np.round(En,4))

Output 1

Shooting Method

The shooting method is a numerical technique to solve boundary value problems by converting them into initial value problems. In this code:

The TISE is transformed into a system of first-order ordinary differential equations (ODEs).
The boundary conditions are specified (i.e., $ \Psi(0) = 0 $ and $ \Psi(L) = 0 $).

The differential equations derived from the TISE are:

Let $ \Psi(x) = y $ (the wave function).

The first derivative $ \frac{d \Psi}{dx} = z $.

Thus, the system becomes:

\[ \frac{dy}{dx} = z \]

\[ \frac{dz}{dx} = \frac{2}{\hbar^2}(V(x) - E) y \]

This can be represented in the code as:

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

u[0] = 0: This sets the initial value of the wave function $ \Psi(x) $ at the starting position (which is usually $ x = 0 $ in a 1D box). Setting $ \Psi(0) = 0 $ means that the wave function starts at zero at this boundary.

u[1] = 0.01: This sets the initial value of the first derivative of the wave function $ \frac{d \Psi}{dx} $ at the starting position $ x = 0 $. A small positive value (0.01) is chosen here, which indicates that the wave function has a slight slope at the beginning. This is necessary to ensure that the wave function can evolve from this point.

The odeint function from SciPy is used to solve the ODEs numerically:

sol = odeint(f, u, x, args=(E, ))

Boundary Value Calculation

The function psiboundval(E) calculates the wave function at the right boundary $ x = L $:

def psiboundval(E):
    sol = odeint(f, u, x, args=(E, ))
    return sol[:, 0][-1]  # Return the wave function at x = L

Visualization

Finally, the code plots $ \Psi(L) $ against the energy values $ E $:

plt.plot(Erange, Y, label="$\Psi(L)$")
plt.axhline(y=0, color='red', linestyle='--', label="Boundary Condition $\Psi(L)=0$")

Eigenvalue Determination

Identifying Zeros

The eigenvalues are identified by locating the energy values $ E $ for which the wave function $ \Psi(L) $ crosses zero. The crossing points correspond to valid eigenvalues of the system.

Bisection Method

The code uses the bisection method to find these eigenvalues, where it checks for changes in the sign of $ \Psi(L) $:

for i in range(len(Y) - 1):
    if np.sign(Y[i]) != np.sign(Y[i + 1]):
        E_eigen = bisect(psiboundval, Erange[i], Erange[i + 1])
        En.append(E_eigen)

This method effectively narrows down the interval in which an eigenvalue exists by looking for points where the function changes from positive to negative (or vice versa).

TISE Solution