In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh_tridiagonal
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh, norm, onenormest
import sympy as sp
from IPython.display import display

Using https://www.latex4technics.com/ for latex typesetting

We start with the Schrödinger's equation:

$$ \left [ \frac{- \hbar^2}{2m} \frac{\partial}{\partial x^2} + V(x,t) \right ] \Psi (x,t) = E \, \Psi (x,t) $$

where in he case of the infinite square well:

$$ \psi(0) = \psi(L) = 0 $$

This can be represented as a eigen equation:

$$ \hat{H} \Psi = E \Psi $$

Where $ \hat{H} = \left [ \frac{- \hbar^2}{2m} \frac{\partial}{\partial x^2} + V(x,t) \right ] $

To make x dimensionless we substitute $ \overline{x} = x / L $ to get

$$ -\frac{1}{2}\frac{d^2\psi}{\overline{x}^2} + \frac{m L^2 V(\overline{x})}{\hbar^2} \psi = \frac{m L^2 E}{ \hbar^2 } \psi $$

In [2]:
# Set up the symbolic variables
x = sp.symbols('x')
n = sp.symbols('n', integer=True, positive=True)
m, hbar, xbar, E, L = sp.symbols(r'm hbar \bar{x} E, L', real=True, positive=True)
Psi = sp.Function('psi')
V = sp.Function('V')

We will set up for a zero potential for the middle part. So from Schrödinger's equation with no potential we have

$$ \frac{- \hbar^2}{2m} \frac{\partial}{\partial x^2} \Psi (x,t) = E \, \Psi (x,t) $$

Ok we want to write x in dimensionless terms so we write $ \bar{x} = x/ L $ then using the second derivative chain rule we can find $ d^2 \psi / dx^2 $ in terms of $ \bar{x} $. To do this we use the second derivative chain rule.

$$ \frac{d^2 \psi}{d x^2} = \frac{d^2\psi}{d\bar{x}^2} \left ( \frac{d\bar{x}}{dx} \right )^2 + \frac{d\psi}{d\bar{x}} \frac{d^2 \bar{x}}{d x^2} $$

Now because $ \displaystyle \frac{d^2 \bar{x}}{d x^2} = 0 $ we can exclude the right hand term. Therefore

$$ \frac{d^2 \psi}{d x^2} = \frac{d^2\psi}{d\bar{x}^2} \left ( \frac{1}{L} \right )^2 $$

Thus our new Schrödinger equation becomes

$$ \frac{- \hbar^2}{2m L^2} \frac{\partial}{\partial \bar{x}^2} \psi (L \bar{x}) = E \, \psi (L \bar{x}) $$

Now define a new $ \displaystyle \psi $ so that $ \psi_{new} (\bar{x}) = \psi_{orig} (L \bar{x}) $ Then we know that the solutions will be the same and related back to the original with just a scaling of the x argument. (ChatGPT claims this is known as Equivalence of differential equations under diffeomorphism).

$$ \frac{\partial}{\partial \bar{x}^2} \psi_{new} (\bar{x}) = -\frac{2 m L^2 E}{\hbar^2} \, \psi_{new} (\bar{x}) $$

In [3]:
eq = sp.Eq(sp.Derivative(Psi(x), x, x), -(2*m*L**2)* E*Psi(x)/hbar**2)
eq
Out[3]:
$\displaystyle \frac{d^{2}}{d x^{2}} \psi{\left(x \right)} = - \frac{2 E L^{2} m \psi{\left(x \right)}}{\hbar^{2}}$
In [4]:
dsoln = sp.dsolve(eq, Psi(x))
dsoln
Out[4]:
$\displaystyle \psi{\left(x \right)} = C_{1} \sin{\left(\frac{\sqrt{2} \sqrt{E} L \sqrt{m} x}{\hbar} \right)} + C_{2} \cos{\left(\frac{\sqrt{2} \sqrt{E} L \sqrt{m} x}{\hbar} \right)}$
In [5]:
expr = dsoln.rhs.args[0]
expr
Out[5]:
$\displaystyle C_{1} \sin{\left(\frac{\sqrt{2} \sqrt{E} L \sqrt{m} x}{\hbar} \right)}$
In [6]:
constants = sorted(
    [s for s in dsoln.free_symbols if s.name.startswith("C")],
    key=lambda s: s.name
)
C1, C2 = constants
constants
Out[6]:
[C1, C2]
In [7]:
# Rather than using sp.symbols and sp.Wild it is possible to match the sin expression directly without
# regard to what is outside it.  The argument of what you pass to the atoms() method determines which
# types are returned
sin_part = next(iter(expr.atoms(sp.sin)))
arg = sin_part.args[0]

k = arg.coeff(x)
k
Out[7]:
$\displaystyle \frac{\sqrt{2} \sqrt{E} L \sqrt{m}}{\hbar}$
In [8]:
soln = dsoln
soln
Out[8]:
$\displaystyle \psi{\left(x \right)} = C_{1} \sin{\left(\frac{\sqrt{2} \sqrt{E} L \sqrt{m} x}{\hbar} \right)} + C_{2} \cos{\left(\frac{\sqrt{2} \sqrt{E} L \sqrt{m} x}{\hbar} \right)}$

This implies that $ C_2 = 0 $ as it must be zero at the boundaries. It must also be zero at the other boundary where x = L. We also ignore the C1 that is used for now

In [9]:
# Because it must be true that ψ = 0 when x = 0 then we can can conclude that it 
# must also be true that C2 = 0
soln.rhs.subs(x,0)
Out[9]:
$\displaystyle C_{2}$
In [10]:
# From that we know we need only deal with the sine part
sin_part
Out[10]:
$\displaystyle \sin{\left(\frac{\sqrt{2} \sqrt{E} L \sqrt{m} x}{\hbar} \right)}$
In [11]:
k = sin_part.args[0].subs(x, L)
k
Out[11]:
$\displaystyle \frac{\sqrt{2} \sqrt{E} L^{2} \sqrt{m}}{\hbar}$

Now to determine the normalisation factor

In [12]:
integral = sp.integrate(sin_part**2,(x,0,L))
normalisation = integral.subs(sp.sin(k), 0)
normalisation
Out[12]:
$\displaystyle \frac{L}{2}$
In [13]:
psifn = sp.sqrt(1/normalisation)*sin_part
psifn
Out[13]:
$\displaystyle \frac{\sqrt{2} \sin{\left(\frac{\sqrt{2} \sqrt{E} L \sqrt{m} x}{\hbar} \right)}}{\sqrt{L}}$

Now we revert back to $ \psi_{old} (x) = \psi_{new} (L x) $

In [14]:
C = sp.Wild("C")
k = sp.Wild("k")

# match x*L , not just x
match = psifn.match(C * sp.sin(k*x*L))

amplitude = match[C]
frequency = match[k]
display(amplitude)
frequency_factor = match[k]
display(frequency_factor)
# Now calculate our original old ψ
psi_old = amplitude*sp.sin(frequency_factor * x)
psi_old
$\displaystyle \frac{\sqrt{2}}{\sqrt{L}}$
$\displaystyle \frac{\sqrt{2} \sqrt{E} \sqrt{m}}{\hbar}$
Out[14]:
$\displaystyle \frac{\sqrt{2} \sin{\left(\frac{\sqrt{2} \sqrt{E} \sqrt{m} x}{\hbar} \right)}}{\sqrt{L}}$
In [15]:
expr = sp.Eq(frequency*L , sp.pi*n)
expr
Out[15]:
$\displaystyle \frac{\sqrt{2} \sqrt{E} L \sqrt{m}}{\hbar} = \pi n$
In [16]:
energy = sp.solve(expr, E)[0]
energy
Out[16]:
$\displaystyle \frac{\pi^{2} \hbar^{2} n^{2}}{2 L^{2} m}$

So $ \displaystyle E_n = \frac{\pi^{2} \hbar^{2} n^{2}}{2 L^{2} m} $

Now the units of $ \hbar $ are $ M L^2 T^{-1} $ or equivalently $ J \cdot s $

Rearranging $ \displaystyle \frac{m L^2 E_n}{\hbar^2} = \frac{\pi^2 n^2}{2} $

$$ \begin{bmatrix}\frac{1}{\Delta x^2}+mL^2V_1/\hbar^2 & -\frac{1}{2 \Delta x^2} & 0 & 0...\\ -\frac{1}{2 \Delta x^2} & \frac{1}{\Delta x^2}+mL^2V_2/\hbar^2 & -\frac{1}{2 \Delta x^2} & 0... \\ ...& ... & ... & -\frac{1}{2 \Delta x^2}\\...0 & 0 & -\frac{1}{2 \Delta x^2} & \frac{1}{\Delta x^2}+mL^2V_{N-1}/\hbar^2 \\ \end{bmatrix} \begin{bmatrix} \psi_1 \\ \psi_2 \\ ... \\ \psi_{N-1} \end{bmatrix} = \frac{mL^2 E}{\hbar^2} \begin{bmatrix} \psi_1 \\ \psi_2 \\ ... \\ \psi_{N-1} \end{bmatrix} $$

$$ \psi_0 = \psi_N = 0$$

$$ \hbar = 1.054571817...×10^{−34} J⋅s $$ However $ \hbar $ is not needed in our evaluation as our potential is either zero or infinity and dividing it by $ \hbar^2 $ is not expected to make any difference on the matrix side but would in the eigenvalue side.

In [17]:
from scipy.constants import hbar as hbar0

Define what $N$ and $dx$ is

In [18]:
# The grid spacing turned out to be tricky and getting it wrong multiplies the
# eigenvalues by a erroneous factor.  I have made the dx formula more explicit so
# its not left to guesswork.
N = 1000
start = 0
end = 1
dx = (end - start)/N
x = np.linspace(start, end, N)

Define potential $mL^2 V$

In [19]:
# This is currently not used as it is more efficient to make the well as wide as
# the matrix.
def mL2V(x):
    # The square well is zero only between 0 and 1
    return np.where(((x<0) | (x>1)), 1e8, 0)
In [20]:
d = 1/dx**2 * np.ones(N) #+ mL2V(x) # main diagnonal
e = -1/(2*dx**2) * np.ones(len(d)-1)  # off diagonal

offsets=[0,-1, 1]
diagonals = [d, e, e]
sparse_matrix = diags(diagonals, offsets, format='csr')

# 1-norm of the matrix
norm_A = norm(sparse_matrix, ord=1)

# Estimated 1-norm of the inverse matrix
norm_inv_A = onenormest(sparse_matrix)

# Condition number estimation
cond_num_est = norm_A * norm_inv_A
print("Estimated Condition Number (1-norm):", cond_num_est)

Estimated Condition Number (1-norm): 4000000000000.0
In [21]:
#w, v = eigh_tridiagonal(d, e)
w,v = eigsh(sparse_matrix, k=10, which='SM')
In [22]:
fig = plt.figure(figsize=(10,5))
ax = fig.gca()
ax.set_xlim(-0.5, 1.5)
ax.set_ylim(-0.5, 5)
plt.plot(x, mL2V(x), lw=3)
# Add a vertical line at x = 2, from y = 0 to y = 12
plt.vlines(x=0, ymin=0, ymax=12, linestyle='-', lw=3)
plt.vlines(x=1, ymin=0, ymax=12, linestyle='-', lw=3)
plt.title('Infinite Potential Well', fontsize=25)
plt.ylabel('$V/V_0$', fontsize=25)
plt.xlabel('$x/L$', fontsize=25)
plt.grid()
#plt.savefig('v3p1.png', dpi=200)
No description has been provided for this image
In [23]:
print (len(x), len(v.T[0]))

1000 1000
In [24]:
fig = plt.figure(figsize=(10,5))
ax = fig.gca()

ax.set_xlim(-0.25, 1.25)
ax.set_ylim(-0.05, 0.05)

line1, = plt.plot(x, v.T[0])
clr = line1.get_color()
plt.hlines(y=0, xmin=-.25, xmax=0, color=clr, linestyle='-', linewidth=4)
plt.hlines(y=0, xmin=1, xmax=1.25, color=clr, linestyle='-', linewidth=4)
plt.vlines(x=0, ymin=0, ymax=12, linestyle='-', lw=4)
plt.vlines(x=1, ymin=0, ymax=12, linestyle='-', lw=4)
plt.plot(x, v.T[0], label='1st eigenstate')
plt.plot(x, v.T[1], label='2nd eigenstate')
plt.plot(x, v.T[2], label='3rd eigenstate')
plt.title('Eigenstates', fontsize=30)
plt.ylabel(r'$ \psi $', fontsize=25)
plt.xlabel('x/L', fontsize=25)
plt.legend()
plt.grid()
#plt.savefig('v3p2.png', dpi=200)
No description has been provided for this image

Now $ \displaystyle \frac{m L^2 E_n}{\hbar^2} = \frac{\pi^2 n^2}{2} $ so we should be able to check the eigenvalues against $ \displaystyle \frac{\pi^2 n^2}{2} $

In [25]:
plt.bar(np.arange(0, 10, 1), w[0:10])
plt.ylabel(r"$ \frac{mL²E}{\hbar^2} $", fontsize=20)
plt.grid()
No description has been provided for this image
In [26]:
from scipy.constants import pi as pi0
errors = []
for i in range(10):
    predicted = pi0**2*(i+1)**2/2
    numerical = w[i]
    error = abs(numerical - predicted)/predicted
    errors.append(error)
    print (f"n={i+1} numerical = {w[i]:.2f} predicted = {pi0**2*(i+1)**2/2:.2f} rel.error = {error:.4f} " )

n=1 numerical = 4.92 predicted = 4.93 rel.error = 0.0020 
n=2 numerical = 19.70 predicted = 19.74 rel.error = 0.0020 
n=3 numerical = 44.32 predicted = 44.41 rel.error = 0.0020 
n=4 numerical = 78.80 predicted = 78.96 rel.error = 0.0020 
n=5 numerical = 123.12 predicted = 123.37 rel.error = 0.0020 
n=6 numerical = 177.29 predicted = 177.65 rel.error = 0.0020 
n=7 numerical = 241.31 predicted = 241.81 rel.error = 0.0020 
n=8 numerical = 315.18 predicted = 315.83 rel.error = 0.0020 
n=9 numerical = 398.89 predicted = 399.72 rel.error = 0.0021 
n=10 numerical = 492.45 predicted = 493.48 rel.error = 0.0021
In [ ]: