%matplotlib widget

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh_tridiagonal
from scipy.integrate import quad
import sympy as sp
import ipywidgets as widgets
import matplotlib.animation as animation

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_2 e^{\alpha x} = C_3 \sin(k x) + C_4 \cos( k x ) $ for $ x = 0 $ and $ C_3 \sin(k x) + C_4 \cos( k x ) = C_1 e^{-\alpha x} $ for x = L:

# 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} $$

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}}$

eq1 = sp.Eq(sp.Derivative(Psi(x), x, x), (alpha)**2*Psi(x))
eq1

$\displaystyle \frac{d^{2}}{d x^{2}} \psi{\left(x \right)} = \alpha^{2} \psi{\left(x \right)}$

# 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}$

right_psi, left_psi = dsoln1.rhs.args
left_psi, right_psi

(C2*exp(\alpha*x), C1*exp(-\alpha*x))

eq2 = sp.Eq(sp.Derivative(Psi(x), x, x), -(k)**2*Psi(x))
eq2

$\displaystyle \frac{d^{2}}{d x^{2}} \psi{\left(x \right)} = - k^{2} \psi{\left(x \right)}$

# schrodinger solutions for center
dsoln = sp.dsolve(eq2, Psi(x))
dsoln2 = dsoln.subs({
    C1: C3,
    C2: C4,
})
dsoln2

$\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

x0 = 0
x1 = L

eqn1 = dsoln2.rhs.subs(x, x0) - left_psi.subs(x,x0)
eqn1

$\displaystyle - C_{2} + C_{4}$

psi_2_deriv = sp.diff(dsoln2.rhs, x)
sp.Eq(sp.diff(dsoln2.lhs,x), psi_2_deriv)

$\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)}$

eqn2 = psi_2_deriv.subs(x, x0) - sp.diff(left_psi,x).subs(x,x0)
sp.Eq(eqn2,0)

$\displaystyle - C_{2} \alpha + C_{3} k = 0$

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

eqn3 = dsoln2.rhs.subs(x,x1) - right_psi.subs(x,x1)
sp.Eq(eqn3, 0)

$\displaystyle - C_{1} e^{- L \alpha} + C_{3} \sin{\left(L k \right)} + C_{4} \cos{\left(L k \right)} = 0$

# 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)

$\displaystyle C_{1} \alpha e^{- L \alpha} + C_{3} k \cos{\left(L k \right)} - C_{4} k \sin{\left(L k \right)} = 0$

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)}$

$\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]$

det = M.det()
det

$\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.

detsimp = sp.cancel(det/(sp.exp(-L*alpha)))
sp.Eq(detsimp, 0)

$\displaystyle \alpha^{2} \sin{\left(L k \right)} + 2 \alpha k \cos{\left(L k \right)} - k^{2} \sin{\left(L k \right)} = 0$

cos_solve = sp.solve(detsimp, sp.cos(L*k))[0]
sp.Eq(sp.cos(L*k) , cos_solve)

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

M0 = M.subs({ sp.cos(L*k): cos_solve})
M0

$\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]$

M0.det()

$\displaystyle 0$

null = M0.nullspace()[0]
null

$\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]$

constants = (null*2*alpha*k*L**2).expand().subs([(L*k, K0), (L*alpha, Alpha0)])
c1, c2, c3, c4 = constants
constants

$\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]$

Calculate the normalisation constant

t1 = c2**2**sp.integrate(sp.exp(Alpha0*x)**2, (x, -sp.oo, x0))
t1

$\displaystyle \left(2 \alpha_{0} k_{0}\right)^{2^{\frac{1}{2 \alpha_{0}}}}$

t2 = c1**2*sp.integrate(sp.exp(-Alpha0*x)**2, (x, x1/L, sp.oo))
t2

$\displaystyle \frac{\left(\alpha_{0}^{2} e^{\alpha_{0}} \sin{\left(k_{0} \right)} + k_{0}^{2} e^{\alpha_{0}} \sin{\left(k_{0} \right)}\right)^{2} e^{- 2 \alpha_{0}}}{2 \alpha_{0}}$

t3 = sp.integrate((c3*sp.sin(K0*x) + c4*sp.cos(K0*x))**2, (x, x0/L, x1/L))
t3

$\displaystyle 2 \alpha_{0}^{4} \sin^{2}{\left(k_{0} \right)} + 2 \alpha_{0}^{4} \cos^{2}{\left(k_{0} \right)} - \frac{2 \alpha_{0}^{4} \sin{\left(k_{0} \right)} \cos{\left(k_{0} \right)}}{k_{0}} + 4 \alpha_{0}^{3} \sin^{2}{\left(k_{0} \right)} + 2 \alpha_{0}^{2} k_{0}^{2} \sin^{2}{\left(k_{0} \right)} + 2 \alpha_{0}^{2} k_{0}^{2} \cos^{2}{\left(k_{0} \right)} + 2 \alpha_{0}^{2} k_{0} \sin{\left(k_{0} \right)} \cos{\left(k_{0} \right)}$

normalisation = sp.sqrt(sp.cancel(t1 + t2 + t3)).subs(L,1)
normalisation

$\displaystyle \frac{\sqrt{2} \sqrt{2 \cdot 2^{2^{\frac{1}{2 \alpha_{0}}}} \alpha_{0} \alpha_{0}^{2^{\frac{1}{2 \alpha_{0}}}} k_{0} k_{0}^{2^{\frac{1}{2 \alpha_{0}}}} + 4 \alpha_{0}^{5} k_{0} \sin^{2}{\left(k_{0} \right)} + 4 \alpha_{0}^{5} k_{0} \cos^{2}{\left(k_{0} \right)} - 4 \alpha_{0}^{5} \sin{\left(k_{0} \right)} \cos{\left(k_{0} \right)} + 9 \alpha_{0}^{4} k_{0} \sin^{2}{\left(k_{0} \right)} + 4 \alpha_{0}^{3} k_{0}^{3} \sin^{2}{\left(k_{0} \right)} + 4 \alpha_{0}^{3} k_{0}^{3} \cos^{2}{\left(k_{0} \right)} + 4 \alpha_{0}^{3} k_{0}^{2} \sin{\left(k_{0} \right)} \cos{\left(k_{0} \right)} + 2 \alpha_{0}^{2} k_{0}^{3} \sin^{2}{\left(k_{0} \right)} + k_{0}^{5} \sin^{2}{\left(k_{0} \right)}}}{2 \sqrt{\alpha_{0}} \sqrt{k_{0}}}$

Find the relationship between $ \alpha $ and k

Find the relationship between $ \alpha_0 = L \alpha $ and $ k_0 = L k $

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_0 $ given $ k_0 $.

expr = sp.expand(detsimp*L**2)
newexpr = expr.subs([(L*k, k0), (L*alpha, alpha0)])
display(newexpr)
f = sp.lambdify([alpha0, k0], newexpr)
normalfn = sp.lambdify( [Alpha0, K0], normalisation)

$\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)}$

from scipy.optimize import root_scalar

xk = []
yalpha = []
normals = []
for lk in np.arange(0.1,12,.1):
    root = (root_scalar(f, args=(lk,), x0=1, method='secant', rtol=1e-3))
    xk.append(lk)
    if root.root > 0:
        yalpha.append(root.root)
        #print(root.root,lk)
        normal = normalisation.subs([(Alpha0, root.root), (K0, lk), (L,1)])
        normals.append(normal.evalf())
    else:
        yalpha.append(np.nan)
        normals.append(np.nan)

%matplotlib inline

from numpy import sin, cos



# 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,15)
plt.plot(xk, y_masked, lw=2)
plt.title(r'$ \alpha_0 $ vs $ k_0$', fontsize=25)
plt.xlabel(r'$k_0$', fontsize=20)
plt.ylabel(r'$\alpha_0$', fontsize=20)
plt.grid()

np.interp(4, xk, np.asarray(normals,dtype=float))

np.float64(15.441617037140869)

x0 = -1
x1 = 2
fig = plt.figure(figsize=(10,5))
ax = fig.gca()
ax.set_xlim(x0, x1)
plt.vlines(x=0, ymin=0, ymax=1, linestyle='-', lw=3, color='grey')
plt.vlines(x=1, ymin=0, ymax=1, linestyle='-', lw=3, color='grey')
plt.hlines(y=1, xmin=x0, xmax=0, linestyle='-', lw=3, color='grey')
plt.hlines(y=1, xmin=1, xmax=x1, linestyle='-', lw=3, color='grey')
plt.hlines(y=0, xmin=0, xmax=1, linestyle='-', lw=3, color='grey')

#kv = 4
for kv, clr in ((1, 'green'),(4,'blue'), (7,'red')):
    av = np.interp(kv, 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
    normalisation_constant = np.interp(kv, xk, np.asarray(normals, dtype=float)) # normalfn(av, kv)
    print (normalisation_constant)
    
    
    #ax.set_ylim(-.25, 1.35)
    label = "k = {}".format(kv)
    xv = np.linspace(x0, 0, 20)
    yv = c2*np.exp(av * xv)/normalisation_constant
    plt.plot(xv, yv, color=clr, label=label)
    
    xv = np.linspace(1, x1, 20)
    yv = c1* np.exp(-av * xv)/normalisation_constant
    plt.plot(xv, yv, color=clr)
    
    xv = np.linspace(0, 1, 40)
    yv = (c3* np.sin(kv * xv) + c4*np.cos(kv*xv))/normalisation_constant
    plt.plot(xv, yv, color=clr)

plt.title('Finite Potential Well', fontsize=25)
plt.ylabel('$V/V_0$', fontsize=25)
plt.xlabel('$x/L$', fontsize=25)
plt.legend()
plt.grid()

1.9239538748257319
15.441617037140869
34.03658070845873

def calc_y(kv, xk, yalpha):
    av = np.interp(kv, 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
    #print (c1,c2,c3,c4)
    #print(av, kv)
    if av < 0 or kv <0 or np.isnan(av) : 
        return None, None

    normalisation_constant = np.interp(kv, xk, np.asarray(normals, dtype=float)) 
    if np.isnan(normalisation_constant):
        return None, None
    #print (normalisation_constant)
    
    
    #ax.set_ylim(-.25, 1.35)
    #label = "k = {}".format(kv)
    xv1 = np.linspace(x0, 0, 20)
    yv1 = c2*np.exp(av * xv1)/normalisation_constant
    #plt.plot(xv, yv, color=clr, label=label)
    
    xv2 = np.linspace(1, x1, 20)
    yv2 = c1* np.exp(-av * xv2)/normalisation_constant
    #plt.plot(xv, yv, color=clr)
    
    xv3 = np.linspace(0, 1, 40)
    yv3 = (c3* np.sin(kv * xv3) + c4*np.cos(kv*xv3))/normalisation_constant
    #plt.plot(xv, yv, color=clr)
    xv = np.concatenate((xv1,xv3, xv2))
    yv = np.concatenate((yv1, yv3, yv2))

    return (xv, yv)

calc_y(1, xk, yalpha)

(array([-1.        , -0.94736842, -0.89473684, -0.84210526, -0.78947368,
        -0.73684211, -0.68421053, -0.63157895, -0.57894737, -0.52631579,
        -0.47368421, -0.42105263, -0.36842105, -0.31578947, -0.26315789,
        -0.21052632, -0.15789474, -0.10526316, -0.05263158,  0.        ,
         0.        ,  0.02564103,  0.05128205,  0.07692308,  0.1025641 ,
         0.12820513,  0.15384615,  0.17948718,  0.20512821,  0.23076923,
         0.25641026,  0.28205128,  0.30769231,  0.33333333,  0.35897436,
         0.38461538,  0.41025641,  0.43589744,  0.46153846,  0.48717949,
         0.51282051,  0.53846154,  0.56410256,  0.58974359,  0.61538462,
         0.64102564,  0.66666667,  0.69230769,  0.71794872,  0.74358974,
         0.76923077,  0.79487179,  0.82051282,  0.84615385,  0.87179487,
         0.8974359 ,  0.92307692,  0.94871795,  0.97435897,  1.        ,
         1.        ,  1.05263158,  1.10526316,  1.15789474,  1.21052632,
         1.26315789,  1.31578947,  1.36842105,  1.42105263,  1.47368421,
         1.52631579,  1.57894737,  1.63157895,  1.68421053,  1.73684211,
         1.78947368,  1.84210526,  1.89473684,  1.94736842,  2.        ]),
 array([0.32886128, 0.33845421, 0.34832698, 0.35848774, 0.36894489,
        0.37970707, 0.39078319, 0.4021824 , 0.41391413, 0.42598808,
        0.43841422, 0.45120284, 0.4643645 , 0.4779101 , 0.49185081,
        0.50619819, 0.52096407, 0.53616068, 0.55180058, 0.5678967 ,
        0.5678967 , 0.57566413, 0.5830531 , 0.59005875, 0.59667649,
        0.60290195, 0.60873105, 0.61415996, 0.6191851 , 0.62380317,
        0.62801114, 0.63180624, 0.63518597, 0.63814811, 0.64069072,
        0.64281212, 0.64451092, 0.645786  , 0.64663653, 0.64706193,
        0.64706195, 0.64663657, 0.64578607, 0.64451101, 0.64281224,
        0.64069087, 0.63814828, 0.63518617, 0.63180646, 0.62801139,
        0.62380345, 0.6191854 , 0.61416029, 0.60873141, 0.60290233,
        0.59667689, 0.59005918, 0.58305354, 0.5756646 , 0.5678972 ,
        0.56789612, 0.55180002, 0.53616013, 0.52096354, 0.50619767,
        0.49185031, 0.47790961, 0.46436403, 0.45120238, 0.43841377,
        0.42598764, 0.41391371, 0.40218199, 0.39078279, 0.37970668,
        0.36894451, 0.35848737, 0.34832662, 0.33845387, 0.32886094]))

%matplotlib widget

plt.close("all")

x0 = -1
x1 = 2
fig = plt.figure(figsize=(10,5))
ax = fig.gca()
ax.set_xlim(x0, x1)
ax.set_ylim(-2,2)
plt.vlines(x=0, ymin=0, ymax=1, linestyle='-', lw=3, color='grey')
plt.vlines(x=1, ymin=0, ymax=1, linestyle='-', lw=3, color='grey')
plt.hlines(y=1, xmin=x0, xmax=0, linestyle='-', lw=3, color='grey')
plt.hlines(y=1, xmin=1, xmax=x1, linestyle='-', lw=3, color='grey')
plt.hlines(y=0, xmin=0, xmax=1, linestyle='-', lw=3, color='grey')

kv_slider = widgets.FloatSlider(
    value=1.0, min=0.3, max=12.0, step=0.05,
    description="kv",
    continuous_update=False
)

def update(change):
    try:
        kv = change["new"]
        #print("change ",change)
        #print("alpha:", av)
        #print("kv:", kv)
        #print("isfinite alpha:", np.isfinite(av))
        #print("isfinite kv:", np.isfinite(kv))
        
        # recalculate y using your finite-square-well formula
    
        xv, yv = calc_y(kv, xk, yalpha)
        try:
            if yv is None:
                line.set_data([], [])
            else:
                line.set_data(xv, yv)
        except ValueError:
            pass
        av = np.interp(kv, xk, yalpha)
        ax.set_title(fr"k = {kv:.3f} $ \alpha $ = {av:.3f}")
        fig.canvas.draw_idle()
    except Exception:
        import traceback
        traceback.print_exc()


kv_initial = 4
xv, yv = calc_y(kv_initial, xk, yalpha)
line, = ax.plot(xv, yv)

kv_slider.observe(update, names="value")
kv_slider.value = kv_initial
display(kv_slider)

plt.title('Finite Potential Well', fontsize=25)
plt.ylabel('$V/V_0$', fontsize=25)
plt.xlabel('$x/L$', fontsize=25)
plt.grid()

FloatSlider(value=4.0, continuous_update=False, description='kv', max=12.0, min=0.3, step=0.05)

Figure

import numpy as np
import matplotlib.animation as animation

kv_values = np.linspace(0.1, 11.5, 1500)
fig, ax = plt.subplots(figsize=(10, 5))
ax.set_xlim(x0, x1)
ax.set_ylim(-2,2)

# draw your well once
ax.vlines(x=0, ymin=0, ymax=1, linestyle='-', lw=3, color='grey')
ax.vlines(x=1, ymin=0, ymax=1, linestyle='-', lw=3, color='grey')
ax.hlines(y=1, xmin=x0, xmax=0, linestyle='-', lw=3, color='grey')
ax.hlines(y=1, xmin=1, xmax=x1, linestyle='-', lw=3, color='grey')
ax.hlines(y=0, xmin=0, xmax=1, linestyle='-', lw=3, color='grey')

xv, yv = calc_y(kv_values[0], xk, yalpha)
#print (xv, yv)
line, = ax.plot(xv, yv)

ax.set_title("Finite Potential Well")
ax.set_ylabel("$V/V_0$", fontsize=25)
ax.set_xlabel("$x/L$", fontsize=25)
ax.grid()

def animate(i):
    if i % 200 == 0:
        print(f"{i}/{len(kv_values)}")
    kv = kv_values[i]
    xv, yv = calc_y(kv, xk, yalpha)

    if yv is not None: 
        #print (max(xv), max(yv))
        line.set_data(xv, yv)
    else:
        line.set_data([], [])
    av = np.interp(kv, xk, yalpha)
    ax.set_title(fr"Finite Potential Well, k = {kv:.3f} $ \alpha $ = {av:.3f}")

    return line,

ani = animation.FuncAnimation(
    fig,
    animate,
    frames=len(kv_values),
    interval=50,
    blit=False
)

ani.save("finite_square_well.mp4", writer="ffmpeg", fps=30, dpi=150)
print("done")

0/1500
0/1500
200/1500
400/1500
600/1500
800/1500
1000/1500
1200/1500
1400/1500
done

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

N = 100
dx = 1/N
x = np.linspace(0, 1, N+1)

Define potential $mL^2 V$

def mL2V(y):
    return np.where(((y<0.25) | (y>0.75)), 500, 0)

d = 1/dx**2 + mL2V(x)[1:-1]
e = -1/(2*dx**2) * np.ones(len(d)-1)

w, v = eigh_tridiagonal(d, e)

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)

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)

plt.bar(np.arange(0, 10, 1), w[0:10])
plt.ylabel('$mL^2 E/ℏ^2$', fontsize=20)
plt.grid()

from sympy import Matrix, exp, sin, cos, integrate, oo, Eq, Rational, lambdify

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

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()

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.

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))

fn(4)

Numerical Solution

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

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 $$

def rootfn(x, r0):
    return x**2/np.cos(x)**2-r0**2

x1 = np.linspace(0, np.pi/2, 10)
y1 = rootfn(x1, 1)
print (y1)

rootval = root_scalar(rootfn, args=(2,), method='bisect', bracket=(0, np.pi/2), rtol=0.0001)
rootval

rootfn(rootval.root, 2)

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()))

        

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()

def fn(x, alpha, k, abfactor):
    if -1 <= x <= 1:
        return abfactor * np.cos(k*x)
    else:
        return np.exp(-alpha * np.abs(x))

def fn2(x, alpha, k, abfactor):
    return fn(x, alpha, k, abfactor)**2

fnv = np.vectorize(fn)

N = 200

xval = np.linspace(-5, 5, N+1)

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)

from IPython.display import display_latex
display_latex("$ 3x $", raw=True)

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 \} $"

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)

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

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()   

# 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)