Navn på gruppens medlemmer:
Numerisk beregning av energier for partikkel i endelig brønnpotensiale¶
For en partikkel i en endelig brønn med bredde $L$ og dybde $V_0$ er $$ V(x) = \begin{cases} -V_0 & \text{for}\,\, 0 < x < L \\ 0 & \text{ellers} \end{cases} $$
Energiene til partikkelen er $$E = \frac{2 \hbar^2 z^2}{m L^2} - V_0,$$ der $z$ er (reelle) løsninger av ligningene $$\tan{(z)} = \sqrt{\frac{z_0^2}{z^2} - 1} \qquad \text{og} \qquad \tan{\left(z + \frac{\pi}{2}\right)} = \sqrt{\frac{z_0^2}{z^2} - 1}$$ med $z_0 = \frac{L}{2 \hbar} \sqrt{2 m V_0}.$ Energiene kan i dette tilfellet ikke angis ved hjelp av enkle analytiske uttrykk. Vi skal her løse disse ligningene numerisk.
# uncomment ONE line to choose matplotlib backend
# if using Jupyter Notebook, use interactive "notebook" backend for best results
# if using Jupyter Lab, use interactive "widget" backend for best results
# if both fail, use static "inline" backend
%matplotlib notebook
#%matplotlib widget
#%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
hbar = 1
Finn energinivåene til en partikkel i en endelig brønn numerisk. Vi har foreslått en framgangsmåte nedenfor.
- Plott funksjonene som involverer $z$ i ligningene over, så du kan se grafisk omtrent hvor løsningene ligger.
- Lokalisér intervaller på $z$-aksen som hver inneholder nøyaktig én løsning.
- Benytt en numerisk metode til å finne løsningen på hvert av intervallene. Kontrollér at de er korrekte ved å plotte dem sammen med funksjonene.
Vi plotter venstre- og høyresidene av ligningene over for $z$ på intervallet $(0, z_0]$ samtidig som vi passer på å unngå singulariteter i $z = 0, \pi/2, \pi, \ldots$. Fra figuren ser vi at annenhver ligning har nøyaktig én løsning på hvert av intervallene mellom $z = 0, \pi/2, \pi, \ldots, z_0$. Vi bruker Newtons metode fra scipy.optimize til å finne løsningene (med startgjetninger på midten av intervallene, der det ikke er noen singulariteter) og bekrefter at de svarer til skjæringspunkter mellom grafene ved å plotte dem sammen med funksjonsuttrykkene.
Vi lager også en figur der vi framstiller energiene som svarer til de numeriske løsningene sammen med brønnpotensialet.
from scipy import optimize
def lhs1(z): return np.tan(z)
def lhs2(z): return np.tan(z + np.pi / 2)
def rhs(z, z0): return np.sqrt((z0/z)**2 - 1)
def f1(z, z0): return lhs1(z) - rhs(z, z0)
def f2(z, z0): return lhs2(z) - rhs(z, z0)
def finwell_energies(L, V0, m):
z0 = L / 2 / hbar * np.sqrt(2 * m * V0)
fig, (ax1, ax2) = plt.subplots(1, 2)
dz = z0 / 1000 # distance between points in plot
# plot well potential, add energies to this plot later
x = np.linspace(-L, +2*L, 1000)
v = np.where(np.logical_and(x > 0, x < L), -V0, 0)
ax2.plot(x, v, color="black", label="V")
ax2.set_xlabel("$x / a_0$")
ax2.set_ylabel("$E / E_h$")
# plot right hand side
z = np.arange(0, z0, dz)[1:] # avoid singularity in z = 0
ax1.plot(z, rhs(z, z0), color="C0", label="$\sqrt{z_0^2/z^2-1}$")
ax1.set_xlabel("$z$")
# loop through all intervals with one solution
zs = []
energies = []
a = 0 # interval left endpoint
b = np.pi/2 # interval right endpoint
i = 0 # interval number
while a < z0:
# equations have solutions on alternating intervals
if i % 2 == 0:
f = f1
lhs = lhs1
lhs_color = "C1"
lhs_label = "$\\tan{(z)}$" if i == 0 else "" # label only once
else:
f = f2
lhs = lhs2
lhs_color = "C2"
lhs_label = "$\\tan{(z + \pi/2)}$" if i == 1 else "" # label only once
# plot separation lines
ax1.axvline(a, color="black", linestyle="dashed")
if b == z0:
# plot last separation line, too
ax1.axvline(b, color="black", linestyle="dashed")
# plot left hand side
z = np.arange(a, b, dz)[1:-1] # slice endpoints to avoid singularities
ax1.plot(z, lhs(z), color=lhs_color, label=lhs_label)
# obtain solution with Newton's method
z_guess = (a + b) / 2
z = optimize.newton(f, z_guess, args=(z0,))
energy = z**2 * hbar**2 / (2 * m * (L/2)**2) - V0
zs.append(z)
energies.append(energy)
# plot solutions
ax1.plot([z], [rhs(z, z0)], "ko")
ax2.plot([0, L], [energy, energy])
ax2.text(+2.1*L/2, energy, f"$E = {energy:.2f} \, E_h$")
# update interval for next iteration
a = a + np.pi/2
b = min(a + np.pi/2, z0)
i = i + 1
# crop meaningful part of plot
ax1.set_ylim(0, 2 * rhs(zs[0], z0))
# make pretty ticks at separation lines
xtick_values = []
xtick_labels = []
for i in range(0, len(zs)):
xtick_values.append(i*np.pi/2)
if i == 0:
xtick_labels.append("0")
elif i % 2 == 1:
xtick_labels.append(f"${i} \pi / 2$")
else: # i % 2 == 0
j = i // 2
xtick_labels.append(f"${j} \pi$")
xtick_values.append(z0)
xtick_labels.append(f"$z_0$")
ax1.set_xticks(xtick_values)
ax1.set_xticklabels(xtick_labels)
ax1.legend(loc="upper right")
ax2.legend(loc="lower left")
m = 1
V0 = 10
L = 3
finwell_energies(L, V0, m)
Shooting-metoden¶
I den første numeriske øvingen løste vi den tidsuavhengige Schrödingerligningen numerisk ved å finne egenverdiene og egenvektorene til en matrise. Vi skal her se på en alternativ og mer illustrerende numerisk teknikk for å løse ligningen. Teknikken kan brukes på alle potensialer, men vi begrenser oss her til potensialer som er symmetriske om $x = 0$, noe som gjør den litt enklere i bruk.
Vi kommer til å få bruk for noen viktige egenskaper til de stasjonære tilstandene.
- Den $n$-te eksiterte tilstanden $\psi(x)$ i et symmetrisk potensial $V(x) = V(-x)$ er
- symmetrisk og krysser ikke $x$-aksen i origo for partalls-$n$. Altså er $\psi(-x) = \psi(x)$ og $\psi(0) \neq 0$.
- antisymmetrisk og krysser $x$-aksen i origo for oddetalls-$n$. Altså er $\psi(-x) = -\psi(x)$ og $\psi(0) = 0$.
- Den $n$-te eksiterte tilstanden $\psi(x)$Â har $n$ nullpunkter.
Schrödingerligningen $$ \hat{H} \psi = -\frac{\hbar^2}{2 m} \psi'' + V \psi = E \psi $$ har, matematisk sett, løsninger for alle verdier av $E$. Gitt et par randverdier og en energi $E$, er det i teorien bare å integrere i vei for å finne en passende funksjon $\psi(x)$. Men i fysikken er vi kun interessert i fysisk akseptable løsninger av ligningen som kan beskrive fysiske partikler. Disse skal beskrive en sannsynlighetstetthet $|\psi|^2$ for posisjonen til en partikkel ved en måling, og må derfor være normerte, dvs. $\int |\psi|^2 \mathrm{d} x$ = 1, noe som krever at $\psi(x) \rightarrow 0$ (tilstrekkelig raskt) når $x \rightarrow \pm \infty$. Vi skal her se hvordan nettopp dette kravet legger (store) begrensninger på hva energiene $E$ kan være.
For å finne fysisk akseptable verdier for $E$ (og dermed $\psi)$, kan vi derfor gå fram på følgende måte:
- Gjett en energi $E$.
- Finn to kjente randverdier for $\psi(x)$.
- Integrér Schrödingerligningen for å finne den passende funksjonen $\psi(x)$.
- Undersøk om $\psi(x)$ er en fysisk akseptabel løsning.
- Hvis $\psi(x) \rightarrow 0$ når $x \rightarrow \pm \infty$, kan $\psi(x)$ normeres til en fysisk akseptabel løsning, og vi har gjettet en riktig verdi for $E$.
- Hvis ikke $\psi(x) \rightarrow 0$ når $x \rightarrow \pm \infty$, er $\psi(x)$ ikke en fysisk akseptabel løsning, og vi har gjettet en gal verdi for $E$.
Kort oppsummert gjetter vi en energi $E$ og "skyter" etter en passende funksjon $\psi(x)$ fra kjente randverdier, før vi undersøker om vi "treffer" $\psi = 0$ i $x = \pm \infty$. Metoden kalles derfor gjerne the shooting method på engelsk.
I den første numeriske øvingen tilnærmet vi den tidsuavhengige Schrödingerligningen som $$-\frac{\hbar^2}{2 m} \frac{\psi_{i+1} - 2 \psi_i + \psi_{i-1}}{{\Delta x}^2} + V_i \psi_i = E \psi_i$$ Vi lar her $i = 0$ svare til $x = 0$. Så lenge vi kjenner to verdier av $\psi_i$ ved siden av hverandre, kan vi bruke denne ligningen til å finne verdiene $\psi_i$ for alle $i$.
I et symmetrisk potensial kan vi for oddetalls-$n$ utnytte $\psi(0) = 0$ til å sette $\psi_0 = 0$ og $\psi_1$ til en vilkårlig verdi forskjellig fra $0$, for eksempel $\psi_1 = 1$. Deretter kan vi benytte den tilnærmede Schrödingerligningen til å finne resten av verdiene $\psi_i$. Hvis det ser ut til at $\psi \rightarrow 0$ for store $|x|$, kan vi normere $\psi$ for å finne en fysisk akseptabel løsning. Den konkrete verdien vi velger for $\psi_1$ er uten betydning, da den uansett blir korrigert ved normeringen.
For partalls-$n$ utnytter vi tilsvarende $\psi(0) \neq 0$ til å sette for eksempel $\psi_0 = 1$. Vi kan så utnytte symmetrien $\psi_1 = \psi_{-1}$ til å finne $\psi_1$ fra den tilnærmede Schrödingerligningen. Så kan vi igjen integrere ligningen og sjekke om energien svarer til en akseptabel løsning.
Skriv en funksjon som skyter etter en symmetrisk eller antisymmetrisk kandidat for energiegenfunksjonen $\psi(x)$ med energi $E$ til en partikkel med masse $m$ i et symmetrisk potensial $V(x)$. Den skal returnere integrasjonspunkter $x_i$, funksjonsverdiene $\psi_i = \psi(x_i)$ i disse punktene og antall nullpunkter i $\psi(x)$. Funksjonen skal skyte ut fra origo i begge retninger og må stoppe å skyte etter en endelig avstand eller når det er åpenbart at $\psi(x)$ divergerer, men den behøver ikke normere $\psi(x)$.
Vi setter opp randverdiene som beskrevet over. For de symmetriske funksjonene finner vi ved å sette $\psi_1 = \psi_{-1}$ i Schrödingerligningen at $$\psi_1 = \psi_0 (1 + m \Delta x^2 (V_0-E) / \hbar^2)$$ Vi løser så ut Schrödingerligningen for $$\psi_{i+1} = \psi_i (2 - 2 m \Delta x^2 (V_i-E) / \hbar^2) - \psi_{i-1}$$ og bruker dette til å integrere den så lenge både funksjonsverdien og avstanden holder seg under en veldig høy verdi. Vi skyter kun ut til den ene siden og teller nullpunkter dobbelt hver gang $\psi_i$ skifter fortegn, og utnytter til slutt (anti)symmetrien til å konstruere hele bølgefunksjonen.
def shoot(V, E, m, is_even, dx=0.001):
nodes = 0
if is_even:
# symmetric wave function
psi0 = 1 # arbitrary
psi1 = psi0 * (1 + m*dx**2/hbar**2*(V(0)-E)) # from Schrödinger equation
else:
# antisymmetric wave function
psi0 = 0
psi1 = 1 # arbitrary
nodes += 1 # node at x = 0
psis = [psi0, psi1]
xs = [0, dx]
x = 2*dx
while np.abs(psi1) < 1e5 and np.abs(x) < 1e5: # probably diverging if wavefunction gets very big
# calculate next value from two previous values
psi1, psi2 = psis[-2], psis[-1]
psi3 = 2 * psi2 - psi1 + psi2 * (V(x-dx) - E) * (2*m*dx**2/hbar**2)
# count node when wave function changes sign
if psi3 * psi2 < 0:
nodes += 2 # one node for both positive and negative x
xs.append(x)
psis.append(psi3)
x += dx
xs = np.array(xs)
psis = np.array(psis)
# use symmetric or antisymmetric property to construct whole wavefunction
xs_rev = xs[1:][::-1]
psis_rev = psis[1:][::-1]
xs = np.concatenate((-xs_rev, xs))
if is_even:
psis = np.concatenate((+psis_rev, psis))
else: # is odd
psis = np.concatenate((-psis_rev, psis))
return xs, psis, nodes
Test funksjonen på et elektron i grunntilstanden i en harmonisk oscillator ved å skyte med energier rett over og under $E = 0.50 \hbar \omega$, for eksempel $E = 0.49 \hbar \omega$ og $E = 0.51 \hbar \omega$. Hva skjer?
Vi ser at $\psi(x)$ divergerer i motsatte retninger for energier rett over og under $E = 0.50 \hbar \omega$. Energien $E = 0.50 \hbar \omega$ er punktet der funksjonen vi skyter bytter fra å divergere mot $+ \infty$ til å divergere mot $- \infty$, og der antall nullpunkter hopper fra 0 til 2. Dette bekrefter at grunntilstandsenergien $E = 0.50 \hbar \omega$ svarer til en fysisk akseptabel løsning av den tidsuavhengige Schrödingerligningen.
m = 1
w = 1
def V(x):
return 1/2*m*w**2*x**2
fig, ax = plt.subplots()
dx = 0.001
for i, E in enumerate(np.linspace(0.49, 0.51, 5)*hbar*w):
x, wave, n = shoot(V, E, m, True)
ax.axhline(0, color="black", linestyle="dashed")
ax.set_xlabel("$x / a_0$")
ax.set_ylabel("$\psi(x)$")
ax.set_ylim(-5, 5)
E_coeff = E / (hbar*w)
ax.plot(x, wave, label=f"$E = {E_coeff:.3f} \hbar \omega$, {n} nodes")
ax.legend(loc="upper center")
Vi kan systematisere gjettingen av $E$ på en måte som lar oss søke etter energien og energiegenfunksjonen til en vilkårlig $n$-te eksiterte tilstand.
Anta at vi vet med sikkerhet at $E$ ligger mellom $E_1$ og $E_2$ (der $E_2 > E_1$), og la disse energiene svare til funksjonene $\psi_1$ og $\psi_2$ med $n_1$ og $n_2$ nullpunkter. Da er $n_1 \leq n$ og $n_2 > n$. Vi kan så skyte ut en ny funksjon $\psi_0$ med energien $E_0 = (E_1 + E_2) / 2$ og $n_0$ nullpunkter. Hvis $n_0 \leq n$, er $E_0$ en bedre nedre skranke for $E$ enn $E_1$. Hvis $n_0 > n$, derimot, er $E_0$ en bedre øvre skranke til $E$ enn $E_2$. Vi kan dermed bytte ut en av skrankene med $E_0$ og gjenta prosessen for å finne et gradvis mindre intervall der energien ligger. Når intervallet er så lite at vi er fornøyd med en hvilken som helst energiverdi på intervallet, stopper vi prosessen.
Men hvordan finner vi to skranker $E_1$ og $E_2$ til å begynne med, hvis vi ikke aner hva den virkelige energien er? En enkel måte er å sjekke antall nullpunkter $n_0$ som svarer til en vilkårlig $E_0$, for eksempel $E_0 = 0$. Denne energien må være enten en nedre eller øvre skranke for $E$. Vi kan så finne den andre skranken ved å sette $E_1 = -1$ (eller $E_2 = +1$), og så senke (eller øke) denne energien gradvis til vi finner den andre skranken. Dette kan gjøres veldig raskt hvis vi øker avstanden mellom $E_1$ og $E_2$ eksponensielt. Så kan vi sette i gang med halveringsmetoden for å finne energien.
Finn energien til den 60. eksiterte tilstanden for et elektron i en harmonisk oscillator.
Du kan gjøre oppgaven på en generell og "automatisk" måte som fungerer for en vilkårlig $n$ og et vilkårlig symmetrisk potensial $V(x)$. Kombinér i så fall metoden for å finne to initielle skranker $E_1$ og $E_2$ med halveringsmetoden.
Du kan også gjøre oppgaven på en "manuell", men like illustrerende måte. Du kan i så fall skyte og plotte funksjoner med varierende energi $E$ og justere energiene etter visuell inspeksjon av plottene inntil du vet med god presisjon hva energien må være.
Vi beskriver her den mest generelle metoden. Koden taler nok best for seg selv, og består i stor grad av å bruke betingelsen for om en energi $E_0$ er en øvre skranke ($n_0 > n$) eller en nedre skranke ($n_0 \leq n$) til å (1) finne to initielle skranker og (2) bruke halveringsmetoden inntil differansen mellom skrankene er så liten at vi er fornøyde med en energi som ligger mellom dem.
Energien til den 60. eksiterte tilstanden ser ut til å være $E = 60.50 \hbar \omega$, som vi forventer fra det analytiske uttrykket.
For fullstendighetens skyld har vi her også normert $\psi(x)$. Siden $\psi(x)$ divergerer for store $|x|$, "kapper" vi av beinene på $\psi(x)$ der de er nærmest $0$ før vi normerer den, slik at $\psi(x)$ etterligner den faktiske energiegenfunksjonen mest mulig.
def normalize(x, wave):
dx = x[1] - x[0]
norm = np.sum(np.abs(wave)**2) * dx
wave /= np.sqrt(norm)
return wave
def is_lower_bound(V, E0, m, n):
_, _, n0 = shoot(V, E0, m, n % 2 == 0)
return n0 <= n
def is_upper_bound(V, E0, m, n):
return not is_lower_bound(V, E0, m, n)
def cut(x, wave):
# cut wavefunction at the furthest point where it is closest to zero
for i in range(1, len(wave)):
if np.abs(wave[i]) > np.abs(wave[i-1]):
break
x = x[i:-i]
wave = wave[i:-i]
return x, wave
def get_stationary_state(n, V, m, tol=1e-5):
dx = 0.001
E0 = 0
if is_lower_bound(V, E0, m, n):
# E0 is lower bound, search for upper bound
E1 = E0
E2 = +1
while not is_upper_bound(V, E2, m, n):
E2 *= 2
else:
# E0 is upper bound, search for lower bound
E2 = E0
E1 = -1
while not is_lower_bound(V, E1, m, n):
E1 *= 2
while (E2 - E1) > tol:
E0 = (E1 + E2) / 2
if is_lower_bound(V, E0, m, n):
E1 = E0 # E0 is better lower bound
else:
E2 = E0 # E0 is better upper bound
E = (E1 + E2) / 2
x, wave, _ = shoot(V, E, m, n % 2 == 0)
x, wave = cut(x, wave) # remove divergent ends
wave = normalize(x, wave)
return x, wave, E
m = 1
w = 1
fig, ax = plt.subplots()
n = 60
x, wave, E = get_stationary_state(n, V, m)
E_coeff = E / (hbar*w)
ax.set_title(f"$n = {n}$, $E = {E_coeff} \hbar \omega$")
ax.set_xlabel("$x / a_0$")
ax.set_ylabel("$|\psi|^2$")
ax.plot(x, np.abs(wave)**2)
ax.set_ylim(0, 0.5)
