In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh_tridiagonal
from scipy.integrate import quad
import sympy as sp

We are looking for some piecewise solutions to 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) $$

If we multiply through by $ -\frac{2m}{\hbar^2} $ and do some rearranging we end up with:

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

We are assuming the solutions to the finite square well look like sinusoids in the well and decaying exponentials outside of it. That means we are looking for solutions of the form $ e^{-\alpha x} $ outside the well and of the form $ sin(k x) $ or $ cos(k x) $, but not a superposition of both due to potential symmetry, within the well.

If $ E \gt V $ then the solutions must be a decaying exponential. If $ E \lt V $ then the solutions will be sinusoids of various frequencies.

We also need to pay attention to the boundary conditions where $ C_1 e^{\alpha x} = C_2 \sin(k x) + C_3 \cos( k x ) $ for $ x = 0 $ and $ C_2 \sin(k x) + C_3 \cos( k x ) = C_4 e^{-\alpha x} $ for x = L:

In [2]:
# Set up the symbolic variables
x = sp.symbols(r'x')
K0, Alpha0 = sp.symbols(r'k_0 \alpha_0', positive=True)
n = sp.symbols('n', integer=True, positive=True)
v0, m, hbar, xbar, E, L, k, alpha = sp.symbols(r'V_0 m hbar \bar{x} E, L k \alpha', real=True, positive=True)
Psi = sp.Function('psi')
Psi_orig = sp.Function(r'\psi_{orig}')
V = sp.Function('V')

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 - V_0) \, \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.

So outside the well the Schrodinger equation becomes

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

and inside the well becomes

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

so

$$ \alpha^2 = \frac{2 m L^2 E }{\hbar^2} \textrm{ and } k^2 = \frac{2 m L^2 (V_0 - E) }{\hbar^2} $$

with

$$ \alpha^2 + k^2 = \frac{2 m L^2 V_0}{\hbar^2} $$

In [3]:
C1, C2, C3, C4, k0, alpha0 = sp.symbols(r"C1 C2 C3 C4 k_0 \alpha_0")

alpha_value = 2 * m*L**2*E/hbar**2
display(sp.Eq(alpha**2, alpha_value))
k2_value = 2*m*L**2*(E - v0)/hbar**2
display(sp.Eq(k**2 , k2_value))
$\displaystyle \alpha^{2} = \frac{2 E L^{2} m}{\hbar^{2}}$
$\displaystyle k^{2} = \frac{2 L^{2} m \left(E - V_{0}\right)}{\hbar^{2}}$
In [4]:
eq1 = sp.Eq(sp.Derivative(Psi(x), x, x), (alpha)**2*Psi(x))
eq1
Out[4]:
$\displaystyle \frac{d^{2}}{d x^{2}} \psi{\left(x \right)} = \alpha^{2} \psi{\left(x \right)}$
In [5]:
# Schrodinger solutions for left and right boundaries
dsoln1 = sp.dsolve(eq1, Psi(x))
display(dsoln1)
$\displaystyle \psi{\left(x \right)} = C_{1} e^{- \alpha x} + C_{2} e^{\alpha x}$
In [6]:
right_psi, left_psi = dsoln1.rhs.args
left_psi, right_psi
Out[6]:
(C2*exp(\alpha*x), C1*exp(-\alpha*x))
In [7]:
eq2 = sp.Eq(sp.Derivative(Psi(x), x, x), -(k)**2*Psi(x))
eq2
Out[7]:
$\displaystyle \frac{d^{2}}{d x^{2}} \psi{\left(x \right)} = - k^{2} \psi{\left(x \right)}$
In [8]:
# schrodinger solutions for center
dsoln = sp.dsolve(eq2, Psi(x))
dsoln2 = dsoln.subs({
    C1: C3,
    C2: C4,
})
dsoln2
Out[8]:
$\displaystyle \psi{\left(x \right)} = C_{3} \sin{\left(k x \right)} + C_{4} \cos{\left(k x \right)}$

Solve for $ \psi (x) $ where left boundary where x = 0¶

In [9]:
x0 = 0
x1 = L
In [10]:
eqn1 = dsoln2.rhs.subs(x, x0) - left_psi.subs(x,x0)
eqn1
Out[10]:
$\displaystyle - C_{2} + C_{4}$
In [11]:
psi_2_deriv = sp.diff(dsoln2.rhs, x)
sp.Eq(sp.diff(dsoln2.lhs,x), psi_2_deriv)
Out[11]:
$\displaystyle \frac{d}{d x} \psi{\left(x \right)} = C_{3} k \cos{\left(k x \right)} - C_{4} k \sin{\left(k x \right)}$
In [12]:
eqn2 = psi_2_deriv.subs(x, x0) - sp.diff(left_psi,x).subs(x,x0)
sp.Eq(eqn2,0)
Out[12]:
$\displaystyle - C_{2} \alpha + C_{3} k = 0$

solve $ \psi (x) $ at boundary conditions at x = L¶

In [13]:
eqn3 = dsoln2.rhs.subs(x,x1) - right_psi.subs(x,x1)
sp.Eq(eqn3, 0)
Out[13]:
$\displaystyle - C_{1} e^{- L \alpha} + C_{3} \sin{\left(L k \right)} + C_{4} \cos{\left(L k \right)} = 0$
In [14]:
# Solve for Psi'(L) 
eqn4 = sp.diff(dsoln2.rhs,x).subs([(x,x1)]) - sp.diff(right_psi,x).subs(x,x1)
sp.Eq(eqn4, 0)
Out[14]:
$\displaystyle C_{1} \alpha e^{- L \alpha} + C_{3} k \cos{\left(L k \right)} - C_{4} k \sin{\left(L k \right)} = 0$
In [15]:
eqns = [eqn1, eqn2, eqn3, eqn4]
mat_data = []
for eqn in eqns:
    display(eqn)
    matrow = [0, 0, 0, 0]
    for term in eqn.as_ordered_terms():
        for idx, C in enumerate([C1, C2, C3, C4]):
            expr  = sp.Wild("expr")
            pattern = C*expr
            result = term.match(pattern)
            if result:
                matrow[idx] = (result[expr])

    mat_data.append(matrow)

M = sp.Matrix(mat_data)
M
$\displaystyle - C_{2} + C_{4}$
$\displaystyle - C_{2} \alpha + C_{3} k$
$\displaystyle - C_{1} e^{- L \alpha} + C_{3} \sin{\left(L k \right)} + C_{4} \cos{\left(L k \right)}$
$\displaystyle C_{1} \alpha e^{- L \alpha} + C_{3} k \cos{\left(L k \right)} - C_{4} k \sin{\left(L k \right)}$
Out[15]:
$\displaystyle \left[\begin{matrix}0 & -1 & 0 & 1\\0 & - \alpha & k & 0\\- e^{- L \alpha} & 0 & \sin{\left(L k \right)} & \cos{\left(L k \right)}\\\alpha e^{- L \alpha} & 0 & k \cos{\left(L k \right)} & - k \sin{\left(L k \right)}\end{matrix}\right]$
In [16]:
det = M.det()
det
Out[16]:
$\displaystyle \alpha^{2} e^{- L \alpha} \sin{\left(L k \right)} + 2 \alpha k e^{- L \alpha} \cos{\left(L k \right)} - k^{2} e^{- L \alpha} \sin{\left(L k \right)}$

Now we take the determinant of this matrix and assume it has less than full rank. (If it did have full rank it would require $ C_1 = C_2 = C_3 = C_4 = 0 $ which is the trivial solution we want to exclude. Besides it would not be normalisable in the slightest.)

Multiplying the determinant through by any constant should have no effect once we have decided that the determinant should be zero.

In [17]:
detsimp = sp.cancel(det/(sp.exp(-L*alpha)))
sp.Eq(detsimp, 0)
Out[17]:
$\displaystyle \alpha^{2} \sin{\left(L k \right)} + 2 \alpha k \cos{\left(L k \right)} - k^{2} \sin{\left(L k \right)} = 0$
In [18]:
cos_solve = sp.solve(detsimp, sp.cos(L*k))[0]
sp.Eq(sp.cos(L*k) , cos_solve)
Out[18]:
$\displaystyle \cos{\left(L k \right)} = \frac{\left(- \alpha^{2} + k^{2}\right) \sin{\left(L k \right)}}{2 \alpha k}$
In [19]:
M0 = M.subs({ sp.cos(L*k): cos_solve})
M0
Out[19]:
$\displaystyle \left[\begin{matrix}0 & -1 & 0 & 1\\0 & - \alpha & k & 0\\- e^{- L \alpha} & 0 & \sin{\left(L k \right)} & \frac{\left(- \alpha^{2} + k^{2}\right) \sin{\left(L k \right)}}{2 \alpha k}\\\alpha e^{- L \alpha} & 0 & \frac{\left(- \alpha^{2} + k^{2}\right) \sin{\left(L k \right)}}{2 \alpha} & - k \sin{\left(L k \right)}\end{matrix}\right]$
In [20]:
M0.det()
Out[20]:
$\displaystyle 0$
In [21]:
null = M0.nullspace()[0]
null
Out[21]:
$\displaystyle \left[\begin{matrix}- \frac{- \alpha^{2} e^{L \alpha} \sin{\left(L k \right)} - k^{2} e^{L \alpha} \sin{\left(L k \right)}}{2 \alpha k}\\1\\\frac{\alpha}{k}\\1\end{matrix}\right]$
In [22]:
constants = (null*2*alpha*k*L**2).expand().subs([(L*k, K0), (L*alpha, Alpha0)])
c1, c2, c3, c4 = constants
constants
Out[22]:
$\displaystyle \left[\begin{matrix}\alpha_{0}^{2} e^{\alpha_{0}} \sin{\left(k_{0} \right)} + k_{0}^{2} e^{\alpha_{0}} \sin{\left(k_{0} \right)}\\2 \alpha_{0} k_{0}\\2 \alpha_{0}^{2}\\2 \alpha_{0} k_{0}\end{matrix}\right]$

Compute the normalisation constant¶

Find the relationship between $ L \alpha $ and $ L k $

In [23]:
from bisect import bisect_left

def nearest_index(arr, value):
    pos = bisect_left(arr, value)

    if pos == 0:
        return 0
    if pos == len(arr):
        return len(arr) - 1

    before = arr[pos - 1]
    after = arr[pos]

    if value - before <= after - value:
        return pos - 1
    else:
        return pos

#x = [1.0, 2.5, 4.0, 7.2, 10.0]

#i = nearest_index(x, 6.8)
#print(i)        # 3
#print(x[i])     # 7.2

Now we take the determinant expression and use a root finding algorithm to find its zero points. This is to find corresponding values of $ \alpha $ given k.

In [47]:
expr = sp.expand(detsimp*L**2)
newexpr = expr.subs([(L*k, k0), (L*alpha, alpha0)])
display(newexpr)
f = sp.lambdify([alpha0, k0], newexpr)
$\displaystyle \alpha_{0}^{2} \sin{\left(k_{0} \right)} + 2 \alpha_{0} k_{0} \cos{\left(k_{0} \right)} - k_{0}^{2} \sin{\left(k_{0} \right)}$
In [51]:
from scipy.optimize import root_scalar
from numpy import sin, cos

xk = []
yalpha = []
for lk in np.arange(0.1,10,.001):
    root = (root_scalar(f, args=(lk,), x0=1, method='secant', rtol=1e-3))
    xk.append(lk)
    yalpha.append(root.root)

# remove the glitches from the plotting to avoid drawing lines between positive infinity
# and negative infinity
threshold = 1
mask = np.zeros_like(yalpha, dtype=bool)
mask[1:] = np.abs(np.diff(yalpha)) > threshold
y_masked = np.ma.masked_array(yalpha, mask)

fig = plt.figure(figsize=(10,5))
ax = fig.gca()
ax.set_ylim(0,10)
plt.plot(xk, y_masked, lw=2)
plt.title(r'$ \alpha L $ vs $ k L$', fontsize=25)
plt.ylabel(r'$L k$', fontsize=20)
plt.xlabel(r'L $\alpha$', fontsize=20)
plt.grid()
No description has been provided for this image
In [50]:
kv = .5
av = np.interp(.5, xk, yalpha)
print (kv, av)
cv = [ float(constants[idx].subs([(Alpha0, av), (K0, kv)]))  for idx in range(0, len(constants)) ]
print(cv)
c1, c2, c3, c4 = cv

0.5 0.1276712269316947
[0.1450571519679163, 0.1276712269316947, 0.03259988437248857, 0.1276712269316947]
In [44]:
fig = plt.figure(figsize=(10,5))
ax = fig.gca()
ax.set_xlim(-.5, 1.5)
ax.set_ylim(-.25, 1.25)

xv = np.linspace(-1.5, 0, 20)
yv = c1*np.exp(av * xv)
plt.plot(xv, yv, color='red')

xv = np.linspace(1, 2.5, 20)
yv = c2* np.exp(-av * xv)
plt.plot(xv, yv, color='red')

xv = np.linspace(0, 1, 40)
yv = c3* np.sin(kv * xv) + c4*np.cos(*xv)
#plt.plot(xv, yv, color='red')

plt.vlines(x=0, ymin=0, ymax=1, linestyle='-', lw=3)
plt.vlines(x=1, ymin=0, ymax=1, linestyle='-', lw=3)
plt.hlines(y=1, xmin=-1, xmax=0, linestyle='-', lw=3)
plt.hlines(y=1, xmin=1, xmax=2, linestyle='-', lw=3)
plt.hlines(y=0, xmin=0, xmax=1, linestyle='-', lw=3)
plt.title('Finite Potential Well', fontsize=25)
plt.ylabel('$V/V_0$', fontsize=25)
plt.xlabel('$x/L$', fontsize=25)
plt.grid()
No description has been provided for this image

Solving Semi-analytically for Solutions¶

We are assuming the solutions to the finite square well look like sinusoids in the well and decaying exponentials outside of it. That means we are looking for solutions of the form $ e^{-\alpha x} $ outside the well and of the form $ sin(k x) $ or $ cos(k x) $, but not a superposition of both due to potential symmetry, within the well.

By setting $ y = \alpha L/2 $ and $ x = k L / 2 $ we arrive at two solutions which can be solved simultaneously:

$ y = x tan (x) $ and $ x^2 + y^2 = R^2_0 $ where $ R^2_0 = \frac{2 m L^2 V_0}{4 \hbar } $

$$ \begin{bmatrix}\frac{1}{\Delta x^2}+mL^2V_1 & -\frac{1}{2 \Delta x^2} & 0 & 0...\\ -\frac{1}{2 \Delta x^2} & \frac{1}{\Delta x^2}+mL^2V_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} \\ \end{bmatrix} \begin{bmatrix} \psi_1 \\ \psi_2 \\ ... \\ \psi_{N-1} \end{bmatrix} = mL^2 E \begin{bmatrix} \psi_1 \\ \psi_2 \\ ... \\ \psi_{N-1} \end{bmatrix} $$

$$ \psi_0 = \psi_N = 0$$

Define what $N$ and $dy$ is

In [19]:
N = 100
dx = 1/N
x = np.linspace(0, 1, N+1)

Define potential $mL^2 V$

In [20]:
def mL2V(y):
    return np.where(((y<0.25) | (y>0.75)), 500, 0)
In [21]:
d = 1/dx**2 + mL2V(x)[1:-1]
e = -1/(2*dx**2) * np.ones(len(d)-1)
In [22]:
w, v = eigh_tridiagonal(d, e)
In [23]:
plt.figure(figsize=(10,5))
plt.plot(x, mL2V(x), lw=3)
plt.title('Potential', fontsize=30)
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 [24]:
plt.figure(figsize=(10,5))
plt.plot(x[1:-1], v.T[0], label="1")
plt.plot(x[1:-1], v.T[1], label="2")
plt.plot(x[1:-1], v.T[2], label="3")
plt.plot(x[1:-1], v.T[3], label="4")
plt.title('Eigenstates', fontsize=20)
plt.ylabel('ψ', fontsize=15)
plt.xlabel('$x/L$', fontsize=15)
plt.grid()
plt.legend()
#plt.savefig('v3p2.png', dpi=200)
Out[24]:
<matplotlib.legend.Legend at 0x7fee2a1eaba0>
No description has been provided for this image
In [25]:
plt.bar(np.arange(0, 10, 1), w[0:10])
plt.ylabel('$mL^2 E/ℏ^2$', fontsize=20)
plt.grid()
No description has been provided for this image
In [26]:
from sympy import Matrix, exp, sin, cos, integrate, oo, Eq, Rational, lambdify
In [27]:
alpha, L, k, E, hbar, V0, V, m = sp.symbols("alpha L k E hbar V_0 V m", real=True, positive=True)
alpha, L, k, E, hbar, V0, V, m
Out[27]:
(alpha, L, k, E, hbar, V_0, V, m)
In [28]:
fig = plt.figure(figsize=(10,6))

plt.title('Finite Square Well', fontsize=20)
plt.ylim(0, 10)


line_color = None
piranges = [(0, np.pi/2), (np.pi/2, 3*np.pi/2), (3*np.pi/2, 5*np.pi/2),  (5*np.pi/2, 7*np.pi/2), ]
for xlow, xhi in piranges:
    x = np.linspace(xlow+0.01, xhi-.01, 200)
    y = x*np.tan(x)
    if not line_color:
        lines = plt.plot(x,y, label='x tan(x)')
        line_color = lines[-1].get_color()
    else:
        plt.plot(x,y, color=line_color)


plt.ylabel('y = kL / 2', fontsize=15)
plt.xlabel('$ x = \\alpha L/2$', fontsize=15)
plt.grid()
#xmin, xmax, ymin, ymax = plt.axis()

# Plot the circles
theta = np.linspace(0, np.pi, 100)
for r in (2,4,6,8,9):
    xc = r * np.sin(theta)
    yc = r * np.cos(theta)
    plt.plot(xc,yc, label="$x^2 + y^2={}$".format(r*r), linestyle='dotted')

# Plot the asymptotes
ymin, ymax = plt.gca().get_ylim()
xvlines = ( [avline[1] for avline in piranges])
plt.vlines(xvlines, ymin, ymax, linestyle='dashed', color='orange', lw=3)
plt.legend()
Out[28]:
<matplotlib.legend.Legend at 0x7fee26232990>
No description has been provided for this image

It was kind of tricky to work out how many roots I would expect to find and whereabouts I would expect to find them. This is kind of important if we want to bracket the roots in a root finding algorithm.

In [29]:
def fn(r):
    n  = int(r/np.pi)+1
    print("Expect to find {} roots".format(n))
    for i in range(0,n):
        print ("Root {} between {} and {}".format(i+1, i*np.pi, (i+0.5)*np.pi))
In [30]:
fn(4)

Expect to find 2 roots
Root 1 between 0.0 and 1.5707963267948966
Root 2 between 3.141592653589793 and 4.71238898038469
In [31]:
alpha, L, k, E, hbar, V0, x, A, m = sp.symbols("alpha L k E hbar V_0 x A m", real=True, positive=True)
display(hbar)
alpha, L, k, E, hbar, V0, x, A, m
$\displaystyle \hbar$
Out[31]:
(alpha, L, k, E, hbar, V_0, x, A, m)

Now we want to set up a joining of the two wavefunctions $ e^{-\alpha x} $ and $ \cos(kx) $ at the point x = L/2 to be continuous and also we want its derivatives joined be continuous also.

Therefore:

$$ e^{-\alpha x} = \cos(kx) $$

and it's derivative is:

$$ -\alpha e^{-\alpha x} = -k \sin (kx) $$

If we substitute one of the join points x = L/2 we get

$$ e^{-\frac{\alpha L}{2}} = \cos \left( \frac{k L}{2} \right) $$

If we agree to set L = 2 then this simplifies to

$$ e^{- \alpha} = \cos \left( k \right) $$

and

$$ -\alpha e^{-\alpha } = -k \sin (k) $$

In [32]:
m1 = sp.Matrix([[exp(-alpha), -cos(k)],[-alpha*exp(-alpha), k*sin(k)]])
m1
Out[32]:
$\displaystyle \left[\begin{matrix}e^{- \alpha} & - \cos{\left(k \right)}\\- \alpha e^{- \alpha} & k \sin{\left(k \right)}\end{matrix}\right]$
In [33]:
m1.det()
Out[33]:
$\displaystyle - \alpha e^{- \alpha} \cos{\left(k \right)} + k e^{- \alpha} \sin{\left(k \right)}$
In [34]:
mdet = sp.simplify(m1.det()/(sp.exp(-alpha)))
sp.Eq(mdet, 0)
Out[34]:
$\displaystyle - \alpha \cos{\left(k \right)} + k \sin{\left(k \right)} = 0$
In [35]:
mdet.subs(alpha, k*sp.tan(k))
Out[35]:
$\displaystyle k \sin{\left(k \right)} - k \cos{\left(k \right)} \tan{\left(k \right)}$
In [36]:
sp.simplify(_)
Out[36]:
$\displaystyle 0$
In [37]:
slice = m1[0,:]*alpha
m1[0,:] = slice
m1
Out[37]:
$\displaystyle \left[\begin{matrix}\alpha e^{- \alpha} & - \alpha \cos{\left(k \right)}\\- \alpha e^{- \alpha} & k \sin{\left(k \right)}\end{matrix}\right]$
In [38]:
row = m1[1,:] + sp.Matrix(slice)
row
Out[38]:
$\displaystyle \left[\begin{matrix}0 & - \alpha \cos{\left(k \right)} + k \sin{\left(k \right)}\end{matrix}\right]$
In [39]:
m1[1,:] = row
m1
Out[39]:
$\displaystyle \left[\begin{matrix}\alpha e^{- \alpha} & - \alpha \cos{\left(k \right)}\\0 & - \alpha \cos{\left(k \right)} + k \sin{\left(k \right)}\end{matrix}\right]$

However $ \displaystyle \alpha \cos{\left(k \right)} - k \sin{\left(k \right )} = 0 $ so that value ought to be replaced by zero in the matrix also.

In [40]:
m1[1,1]=0
m1
Out[40]:
$\displaystyle \left[\begin{matrix}\alpha e^{- \alpha} & - \alpha \cos{\left(k \right)}\\0 & 0\end{matrix}\right]$

This means $ \begin{pmatrix}e^{-\alpha} & -\cos(k)\\ 0 & 0 \end{pmatrix} \begin{pmatrix}A\\B\end{pmatrix} = \begin{pmatrix}0\\0\end{pmatrix} $ or that $ A e^{- \alpha} -B \cos(k) = 0 $

Therefore we can write $ A = B \cos(k) e^{ \alpha} $ or $ B = A \sec(k) e^{ -\alpha } $ and we can set $ A = 1 $ for the pre normalised value.

Now if assume by the same symmetry argument that half of the wave function probability density function integral must be 0.5 then we should be able to compute both A and B.

In [41]:
m2 = sp.Matrix([[exp(-alpha), -sin(k)],[-alpha*exp(-alpha), -k*cos(k)]])
m2
Out[41]:
$\displaystyle \left[\begin{matrix}e^{- \alpha} & - \sin{\left(k \right)}\\- \alpha e^{- \alpha} & - k \cos{\left(k \right)}\end{matrix}\right]$
In [42]:
m2.det()
Out[42]:
$\displaystyle - \alpha e^{- \alpha} \sin{\left(k \right)} - k e^{- \alpha} \cos{\left(k \right)}$
In [43]:
mdet = sp.simplify(m2.det()/(sp.exp(-alpha)))
sp.Eq(mdet, 0)
Out[43]:
$\displaystyle - \alpha \sin{\left(k \right)} - k \cos{\left(k \right)} = 0$

Compute the normalisation constant for cos(kx)¶

In [44]:
f1 = exp(- alpha * x)
f1
Out[44]:
$\displaystyle e^{- \alpha x}$
In [45]:
f2 = A * cos(k*x)
f2
Out[45]:
$\displaystyle A \cos{\left(k x \right)}$
In [46]:
f3 = A*sin(k*x)
f3
Out[46]:
$\displaystyle A \sin{\left(k x \right)}$
In [47]:
nf1 = integrate(f2**2, (x, 0, 1)) + integrate(f1**2, (x, 1, oo))
nf1
Out[47]:
$\displaystyle \frac{A^{2} \left(\frac{k}{2} + \frac{\sin{\left(k \right)} \cos{\left(k \right)}}{2}\right)}{k} + \frac{e^{- 2 \alpha}}{2 \alpha}$
In [48]:
nf2 = integrate(f3**2, (x, 0, 1)) + integrate(f1**2, (x, 1, oo))
nf2
Out[48]:
$\displaystyle \frac{A^{2} \left(\frac{k}{2} - \frac{\sin{\left(k \right)} \cos{\left(k \right)}}{2}\right)}{k} + \frac{e^{- 2 \alpha}}{2 \alpha}$
In [49]:
e1 = Eq(nf1, Rational(1,2))
e1
Out[49]:
$\displaystyle \frac{A^{2} \left(\frac{k}{2} + \frac{\sin{\left(k \right)} \cos{\left(k \right)}}{2}\right)}{k} + \frac{e^{- 2 \alpha}}{2 \alpha} = \frac{1}{2}$
In [50]:
e2 = Eq(nf2, Rational(1,2))
e2
Out[50]:
$\displaystyle \frac{A^{2} \left(\frac{k}{2} - \frac{\sin{\left(k \right)} \cos{\left(k \right)}}{2}\right)}{k} + \frac{e^{- 2 \alpha}}{2 \alpha} = \frac{1}{2}$
In [51]:
# Compute the normalisation value
normalisation1 = sp.simplify(sp.solveset(e1, A).args[0])
Eq(A, normalisation1)
Out[51]:
$\displaystyle A = \frac{\sqrt{2} \sqrt{k} \sqrt{\frac{\alpha e^{2 \alpha} - 1}{2 k + \sin{\left(2 k \right)}}} e^{- \alpha}}{\sqrt{\alpha}}$
In [52]:
# Compute the normalisation value
normalisation2 = sp.simplify(sp.solveset(e2, A).args[0])
Eq(A, normalisation2)
Out[52]:
$\displaystyle A = \frac{\sqrt{2} \sqrt{k} \sqrt{\frac{\alpha e^{2 \alpha} - 1}{2 k - \sin{\left(2 k \right)}}} e^{- \alpha}}{\sqrt{\alpha}}$

Numerical Solution¶

Now to find the solutions for $ \alpha $ and k which we are going to have to solve for numerically.

In [53]:
from scipy.optimize import root_scalar

$$ \alpha^2 + k^2 = R_0^2 $$

and since we know

$$ \alpha = k \tan k $$

then we can substitute that into our first equation to get:

$$ \left( k^2 \tan^2(k) \right) + k^2 = R_0^2 \Rightarrow k^2( 1 + \tan^2(k)) = k^2 \sec^2(k) = R_0^2 $$

In [54]:
def rootfn(x, r0):
    return x**2/np.cos(x)**2-r0**2
In [55]:
x1 = np.linspace(0, np.pi/2, 10)
y1 = rootfn(x1, 1)
print (y1)

[-1.00000000e+00 -9.68591166e-01 -8.62011436e-01 -6.34459096e-01
 -1.69448083e-01  8.43146052e-01  3.38649084e+00  1.17599051e+01
  6.36538288e+01  6.58079015e+32]
In [56]:
rootval = root_scalar(rootfn, args=(2,), method='bisect', bracket=(0, np.pi/2), rtol=0.0001)
rootval
Out[56]:
      converged: True
           flag: converged
 function_calls: 16
     iterations: 14
           root: 1.0297804776674795
         method: bisect
In [57]:
rootfn(rootval.root, 2)
Out[57]:
np.float64(-0.001813983187332724)
In [58]:
allroots=[]
for r in range(1,10):
    # calculate the number of roots we expect to find for this radius
    n = int(r/np.pi)+1
    roots = []
    for ii in range(n):
        lowrange = ii*np.pi
        hirange = (ii+.5)*np.pi
        try:
            root = root_scalar(rootfn, args=(r,), method='bisect', bracket=(lowrange+0.01, hirange-0.01), x0=lowrange + 0.01, x1=hirange-0.01, rtol=0.0001)
            if not root.converged:
                print("Warning root {} for radius {} did not converge".format(ii+1, r))
            roots.append(root.root)
        except ValueError:
            print (roots)
            print(r, lowrange, hirange)
            raise
    print (r, roots)
    allroots.append((r, roots.copy()))

        

1 [0.7390654728477594]
2 [1.0298871106698615]
3 [1.170162999116519]
4 [1.252321886546815, 3.59513706199805]
5 [1.3064634575262268, 3.837259681867517]
6 [1.344703308357839, 3.9856757365803097]
7 [1.373288543385431, 4.088658305156534, 6.616141671133699]
8 [1.3954373678770082, 4.165138006819758, 6.830436280744518]
9 [1.413042843754929, 4.224201538797298, 6.96825118869211]
In [59]:
fig = plt.figure(figsize=(10,6))

plt.title('Finite Square Well With Discrete Energies', fontsize=15)
plt.ylim(0, 10)


line_color = None
piranges = [(0, np.pi/2), (np.pi/2, 3*np.pi/2), (3*np.pi/2, 5*np.pi/2),  (5*np.pi/2, 7*np.pi/2), ]
for xlow, xhi in piranges:
    x1 = np.linspace(xlow+0.01, xhi-.01, 200)
    y1 = x1*np.tan(x1)
    if not line_color:
        lines = plt.plot(x1,y1, label='k tan(k)')
        line_color = lines[-1].get_color()
    else:
        plt.plot(x1,y1, color=line_color)


plt.xlabel('k', fontsize=15)
plt.ylabel(r'$ \alpha $', fontsize=15)
plt.grid()
#xmin, xmax, ymin, ymax = plt.axis()

# Plot the circles
theta = np.linspace(0, np.pi, 100)
for r in range(1,10):
    xc = r * np.sin(theta)
    yc = r * np.cos(theta)
    plt.plot(xc,yc, color='lightgreen', linestyle='dotted', lw=3)
    
# Plot the intersections of the circles and the x tan x lines.
scatterx = []
scattery = []
energies = []
for r, roots in allroots:
    for root in roots:
        k1 = root
        alpha1 = k1*np.tan(k1)
        scatterx.append(k1)
        scattery.append(alpha1)
        
        abfactor = np.exp(-alpha1)/np.cos(k1)
        print ("r={}, k={:.4f}, α={:.4f}, AB factor={:7.4f} E={:.4f} hbar^2/2m"\
               .format(r, k1,alpha1, abfactor, k1**2))
        energies.append((r, k1,alpha1, abfactor))
plt.scatter(scatterx, scattery, c='red')

# Plot the asymptotes
ymin, ymax = plt.gca().get_ylim()
xvlines = ( [avline[1] for avline in piranges])
plt.vlines(xvlines, ymin, ymax, linestyle='dashed', color='orange', lw=3)
plt.legend()

r=1, k=0.7391, α=0.6736, AB factor= 0.6899 E=0.5462 hbar^2/2m
r=2, k=1.0299, α=1.7146, AB factor= 0.3497 E=1.0607 hbar^2/2m
r=3, k=1.1702, α=2.7628, AB factor= 0.1618 E=1.3693 hbar^2/2m
r=4, k=1.2523, α=3.7984, AB factor= 0.0716 E=1.5683 hbar^2/2m
r=4, k=3.5951, α=1.7524, AB factor=-0.1929 E=12.9250 hbar^2/2m
r=5, k=1.3065, α=4.8268, AB factor= 0.0307 E=1.7068 hbar^2/2m
r=5, k=3.8373, α=3.2038, AB factor=-0.0529 E=14.7246 hbar^2/2m
r=6, k=1.3447, α=5.8459, AB factor= 0.0129 E=1.8082 hbar^2/2m
r=6, k=3.9857, α=4.4832, AB factor=-0.0170 E=15.8856 hbar^2/2m
r=7, k=1.3733, α=6.8624, AB factor= 0.0053 E=1.8859 hbar^2/2m
r=7, k=4.0887, α=5.6822, AB factor=-0.0058 E=16.7171 hbar^2/2m
r=7, k=6.6161, α=2.2881, AB factor= 0.1074 E=43.7733 hbar^2/2m
r=8, k=1.3954, α=7.8759, AB factor= 0.0022 E=1.9472 hbar^2/2m
r=8, k=4.1651, α=6.8356, AB factor=-0.0021 E=17.3484 hbar^2/2m
r=8, k=6.8304, α=4.1620, AB factor= 0.0182 E=46.6549 hbar^2/2m
r=9, k=1.4130, α=8.8829, AB factor= 0.0009 E=1.9967 hbar^2/2m
r=9, k=4.2242, α=7.9543, AB factor=-0.0007 E=17.8439 hbar^2/2m
r=9, k=6.9683, α=5.6936, AB factor= 0.0043 E=48.5565 hbar^2/2m
Out[59]:
<matplotlib.legend.Legend at 0x7fee27ede990>
No description has been provided for this image
In [60]:
def fn(x, alpha, k, abfactor):
    if -1 <= x <= 1:
        return abfactor * np.cos(k*x)
    else:
        return np.exp(-alpha * np.abs(x))
In [61]:
def fn2(x, alpha, k, abfactor):
    return fn(x, alpha, k, abfactor)**2
In [62]:
fnv = np.vectorize(fn)
In [63]:
N = 200

xval = np.linspace(-5, 5, N+1)
In [71]:
plt.figure(figsize=(10,5))
r, k, alp, abfactor = energies[0]

# Normalise the area under the curve
area = quad(fn2, -np.inf, np.inf, args=(alp, k, abfactor) )
yval = fnv(xval, alp, k, abfactor)/np.sqrt(area[0])
#display(exp(-alpha))
#print (area)


plt.plot(xval,yval)
plt.title('Wavefunction', fontsize=20)
plt.ylabel(r'$\Psi$', fontsize=15)
plt.xlabel('$x$', fontsize=15)
plt.grid()

ymin, ymax = plt.gca().get_ylim()
xvlines = [-1, 1]
plt.vlines(xvlines, ymin, ymax, linestyle='dashed', color='orange', lw=3)
#plt.legend()

#plt.savefig('v3p1.png', dpi=200)
Out[71]:
<matplotlib.collections.LineCollection at 0x7fee25da1f90>
No description has been provided for this image
In [72]:
from IPython.display import display_latex
display_latex("$ 3x $", raw=True)
$ 3x $
In [73]:
plt.figure(figsize=(10,5))

#s = r"""$ \begin{{ matrix }}{:.4f} e^{{- {:.4f} x}} & \text{{if }} x \leq 1 \\ {:.4f} \sin({:.4f} x) & \text{{if }} x \ge 1  \end{{matrix}} $
s = r"$ v_0 = %d \frac{\hbar^2}{2m} \Rightarrow \Psi = \left \{ \begin{matrix} %.3f \, e^{- %.3f \, x} & \text{if } |x| \lt 1 " + \
    r"\\ %.3f \cos (%.3f x) & \text{if } |x| \ge 1   \end{matrix} \right \} $"

<Figure size 1000x500 with 0 Axes>
In [74]:
plt.figure(figsize=(10,5))

#s = r"""$ \begin{{ matrix }}{:.4f} e^{{- {:.4f} x}} & \text{{if }} x \leq 1 \\ {:.4f} \sin({:.4f} x) & \text{{if }} x \ge 1  \end{{matrix}} $
#s = r"$ v_0 = %d \frac{\hbar^2}{2m} \Rightarrow \Psi = \left \{ \begin{matrix} %.3f \, e^{- %.3f \, x} & \text{if } |x| \lt 1 " + \
#    r"\\ %.3f \cos (%.3f x) & \text{if } |x| \ge 1   \end{matrix} \right \} $"
#print (repr(s))
eqns = []

for ii in range(0,3):
    r, k, alp, abfactor = energies[ii]
    area = quad(fn2, -np.inf, np.inf, args=(alp, k, abfactor) )
    #print(r, k, alp)
    #print (area)
    s1 = s % (ii+1, 1/area[0], alp, abfactor/area[0], k)
    #print (alp, np.exp(-alp), abfactor*np.cos(k))
    eqns.append((exp(-alp*x)/np.sqrt(area[0]), abfactor*cos(k*x)/np.sqrt(area[0])))
    #print (fnv(1, alp, k , abfactor))
    display_latex(s1, raw=True)
    #display(abfactor/area[0]*sin(k*x))
    
    
    y = fnv(xval, alp, k, abfactor)/np.sqrt(area[0])
    #y = fnv(x, alp, k, abfactor)
    plt.plot(xval,y, label=r"$V_0 ={} \frac{{\hbar^2}}{{2m}}$".format(ii+1))
    #print (r)
plt.title('Finite Square Well - 1st Energy Level - Various Potentials', fontsize=15)
plt.ylabel(r'$\Psi$', fontsize=15)
plt.xlabel('$x$', fontsize=15)
plt.grid()

ymin, ymax = plt.gca().get_ylim()
xvlines = [-1, 1]
plt.vlines(xvlines, ymin, ymax, linestyle='dashed', color='orange', lw=3)
plt.legend()

#plt.savefig('v3p1.png', dpi=200)
$ v_0 = 1 \frac{\hbar^2}{2m} \Rightarrow \Psi = \left \{ \begin{matrix} 0.846 \, e^{- 0.674 \, x} & \text{if } |x| \lt 1 \\ 0.583 \cos (0.739 x) & \text{if } |x| \ge 1 \end{matrix} \right \} $
$ v_0 = 2 \frac{\hbar^2}{2m} \Rightarrow \Psi = \left \{ \begin{matrix} 5.166 \, e^{- 1.715 \, x} & \text{if } |x| \lt 1 \\ 1.806 \cos (1.030 x) & \text{if } |x| \ge 1 \end{matrix} \right \} $
$ v_0 = 3 \frac{\hbar^2}{2m} \Rightarrow \Psi = \left \{ \begin{matrix} 28.036 \, e^{- 2.763 \, x} & \text{if } |x| \lt 1 \\ 4.537 \cos (1.170 x) & \text{if } |x| \ge 1 \end{matrix} \right \} $
Out[74]:
<matplotlib.legend.Legend at 0x7fee25c3c690>
No description has been provided for this image

Now recall:

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

Perform a sanity test on the derived equations

In [75]:
fig = plt.figure(figsize=(10,6))

plt.title('Finite Square Well', fontsize=15)

colors = ['blue','orange','green']

for clr, elmt in zip(colors, eqns):
    expx, cosx = elmt
    
    xseg1 = np.linspace(-5,-1, N)
    f1 = lambdify(x, expx)
    yval = f1(-xseg1)
    plt.plot(xseg1, yval, color=clr)
    
    xseg2 = np.linspace(1,5, N)
    yval = f1(xseg2)
    plt.plot(xseg2, yval, color=clr)
    
    f2 = lambdify(x, cosx)
    xseg3 = np.linspace(-1,1, N)
    yval = f2(xseg3)
    plt.plot(xseg3, yval, color=clr)
    
ymin, ymax = plt.gca().get_ylim()
xvlines = [-1, 1]
plt.vlines(xvlines, ymin, ymax, linestyle='dashed', color='orange', lw=2)
    
plt.grid()
No description has been provided for this image
In [76]:
# The code under here is doubtful given the way it tries to calculate the energies and ends up
# with contradictory answers.
s = ""
for elmt in eqns:
    expx, cosx = elmt
    expxdiff = expx.diff(x,x)
    cosxdiff = cosx.diff(x,x)
    expr_exp = sp.solveset(sp.simplify((expx.diff(x,x)+2*m/hbar**2*(E-V0)*expx)/expx), E)
    cosx_exp = sp.solveset(sp.simplify((cosx.diff(x,x)+2*m/hbar**2*(E-V0)*cosx)/cosx), E)
    #display(Eq(E,expr_exp.args[0]), Eq(E, cosx_exp.args[0]))
    e1latex = (sp.latex(Eq(E,expr_exp.args[0])))
    e2latex = (sp.latex(Eq(E,cosx_exp.args[0])))
    s += "{} & {} ".format(e1latex, e2latex) 
    
    # use this to detect if at last index of for loop
    if elmt is not eqns[-1]:
        s += r" \\ "
    #display_latex ("$ "+e1latex+" $", raw=True)

    #display_latex('E = V_{0} - \\frac{0.226846595732537 \\hbar^{2}}{m}', raw=True)
s = r"$  \left \{ \begin{matrix} " + s + r"\end{matrix} \right \} $"

display_latex(s, raw=True)
$ \left \{ \begin{matrix} E = V_{0} - \frac{0.226846595732537 \hbar^{2}}{m} & E = V_{0} + \frac{0.273108886577841 \hbar^{2}}{m} \\ E = V_{0} - \frac{1.46988327214104 \hbar^{2}}{m} & E = V_{0} + \frac{0.530333730361958 \hbar^{2}}{m} \\ E = V_{0} - \frac{3.81657632689217 \hbar^{2}}{m} & E = V_{0} + \frac{0.684640722250683 \hbar^{2}}{m} \end{matrix} \right \} $
In [ ]: