Navn på gruppens medlemmer:
Tidsutvikling av bølgefunksjonen¶
Utviklingen av en tilstand $\Psi = \Psi(x, t)$ over tid er beskrevet av den tidsavhengige Schrödingerlikningen $$ i \hbar \frac{\partial \Psi}{\partial t} = \hat{H} \Psi = -\frac{\hbar^2}{2m} \frac{\partial^2 \Psi}{\partial x^2} + V \Psi $$
Ved hjelp av separasjon av variable er det mulig å vise at en generell løsning av denne ligningen er $$ \Psi(x, t) = \sum_n c_n \psi_n(x) e^{-i E_n t / \hbar} $$ der summen går over alle de stasjonære tilstandene $\psi_n$ med energier $E_n$ som løser den tidsuavhengige Schrödingerligningen $$ \hat{H} \psi_n = -\frac{\hbar^2}{2 m} \psi_n'' + V \psi_n = E_n \psi_n $$
Det viser seg at enhver tilstand $\Psi(x, t)$ kan skrives som en slik lineærkombinasjon av de stasjonære tilstandene, og vi sier derfor at de utgjør et fullstendig sett. Gitt en vilkårlig starttilstand $\Psi_0 = \Psi(x, 0)$, er det bare snakk om å finne de rette konstantene $c_n$ som oppfyller $$ \Psi_0 = \sum_n c_n \psi_n $$ og vi har da all informasjonen vi trenger for å beregne $\Psi$ for alle $t$. Dette er ingen vanskelig oppgave, da de stasjonære tilstandene er ortogonale, og ortonormale, så lenge vi velger dem normerte, dvs. $$\int \psi_m \psi_n \mathrm{d}x = \delta_{m n} = \begin{cases}0 & \text{for}\,\, m \neq n \\ 1 & \text{for}\,\, m = n\end{cases}$$ Vi kan utnytte dette til å "plukke ut" en av konstantene $c_m$ i lineærkombinasjonen ved å integrere begge sider av likningen: $$ \int \Psi_0 \psi_m \mathrm{d}x = \int \sum_n c_n \psi_n \psi_m \mathrm{d}x = \sum_n c_n \int \psi_n \psi_m \mathrm{d}x = \sum_n c_n \delta_{m n} = c_m $$
I denne øvingen skal vi se nærmere på egenskapene om ortogonalitet og fullstendighet i en diskretisert modell og utnytte dette til å studere tidsutviklingen til vilkårlige starttilstander i vilkårlige potensialer.
Vi bruker samme modell for diskretiseringen av rommet som i den første numeriske øvingen. Der diskretiserte vi rommet i $N + 2$ ekvidistante punkter fra $x_0$ til $x_{N+1}$ og satte $V(x \leq x_0) = V(x \geq x_{N+1}) = \infty$, slik at $\psi_n(x \leq x_0) = \psi_n(x \geq x_{N+1}) = 0$. Vi innførte så et vilkårlig potensial $\boldsymbol{V} = [V(x_1), \ldots, V(x_N)]^T$ på rutenettet $\boldsymbol{x} = [x_1, \ldots, x_N]^T$ og fant de stasjonære tilstandene $\boldsymbol{\psi_n} = [\psi_n(x_1), \ldots, \psi_n(x_N)]^T$ og energiene $E_n$ ved å finne egenverdiene og egenvektorene til Hamiltonmatrisen $H$.
Hamiltonmatrisen $H$ er en reell og symmetrisk matrise, og har da ifølge et teorem fra lineæralgebraen ortogonale egenvektorer. Vi kan så oversette formen for den generelle løsningen på den tidsavhengige Schrödingerlikningen til vår diskretiserte modell med vektorlikningen $$ \boldsymbol{\Psi} = \boldsymbol{\Psi}(x, t) = [\Psi(x_1, t), \ldots, \Psi(x_N, t)]^T = \sum_{n=1}^{N} c_n \boldsymbol{\psi_n} e^{-i E_n t / \hbar} $$ I denne diskretiserte modellen får vi dermed helt tilsvarende ortogonalitets- og fullstendighetsegenskaper som med de eksakte funksjonene $\psi$ og $\Psi$.
# 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 matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np
from IPython.display import HTML, clear_output
# helpful to raise exception upon complex number calculation error
# instead of just a useless warning
import warnings
warnings.filterwarnings(action="error", category=np.ComplexWarning)
For å finne de stasjonære tilstandene $\boldsymbol{\psi_n}$ og energiene $E_n$ til en partikkel i et vilkårlig potensial $\boldsymbol{V}$, bruker vi teknikken fra den første numeriske øvingen.
Vi har lagt ved fungerende funksjonalitet for dette i funksjonen get_stationary_states() under.
Merk at
- det brukes atomære enheter med $\hbar = m_e = 1$, slik at naturlige avstander og tider også har tallverdier i størrelsesorden rundt 1,
- de stasjonære tilstandene $\boldsymbol{\psi_n}$ normeres i forstanden $\boldsymbol{\psi_m} \cdot \boldsymbol{\psi_n} = \sum_{i=1}^N \psi_{m,i} \psi_{n,i} \Delta x = \delta_{mn}$ (der $\psi_{n,i}$ er $i$-te element i vektoren $\boldsymbol{\psi_n}$), i tråd med normeringsintegralet $\int \psi_m \psi_n \mathrm{d}x \approx \sum_{i=1}^N \psi_m(x_i) \psi_n(x_i) \Delta x$, og at
- de stasjonære tilstandene kommer ut i en $N \times N$ matrise med $\boldsymbol{\psi_n}$ langs den $n$-te raden, som så kan hentes ut med
waves[n-1].
Du kan selvfølgelig erstatte dette med din egen løsning fra den første numeriske øvingen, da du kanskje er bedre kjent med denne. For å passe inn med noe inkludert programkode senere i øvingen, anbefaler vi likevel at kravene over er tilfredsstilt.
hbar = 1
def hamiltonian(x, v, m):
dx = x[1] - x[0]
diag = hbar**2 / (m*dx**2) + v
semidiag = np.full(len(x)-1, -hbar**2 / (2*m*dx**2))
H = np.diag(diag, k=0) + np.diag(semidiag, k=+1) + np.diag(semidiag, k=-1)
return H
def normalize(x, wave):
dx = x[1] - x[0]
norm = np.sum(np.abs(wave)**2) * dx
wave /= norm**(1/2)
return wave
def get_stationary_states(x, v, m):
H = hamiltonian(x, v, m)
eigvals, eigvecs = np.linalg.eigh(H) # eigh converges more often than eigh_tridiagonal
energies, waves = eigvals, eigvecs.T
for i in range(0, len(waves)):
waves[i] = normalize(x, waves[i])
return energies, waves
# example usage:
#x = np.linspace(-10, +10, 1000)
#m = 1
#w = 1
#v = 1/2*m*w**2*x**2
#energies, waves = get_stationary_states(x, v, m)
#fig, ax = plt.subplots()
#ax.plot(x, np.abs(waves[0])**2, label=f"$E_0 = {energies[0]:.2f} \, E_h$")
#ax.plot(x, np.abs(waves[1])**2, label=f"$E_1 = {energies[1]:.2f} \, E_h$")
#ax.legend()
Nedenfor har vi konstruert noen potensialer $\boldsymbol{V}$.
For hvert av disse potensialene, finn de stasjonære tilstandene for et elektron og bekreft at de er ortonormale og utgjør et fullstendig sett.
Ta deg friheten til å endre litt på potensialene, om de ikke "passer inn" i programmet ditt i formen under.
Vi beregner $\boldsymbol{\psi_m} \cdot \boldsymbol{\psi_n} = \sum_i \psi_m(x_i) \psi_n(x_i) \Delta x$ for alle par av $m$ og $n$, lagrer verdiene i en matrise og viser denne grafisk. Matrisen ser ut til å være identitetsmatrisen. Det bekrefter at de $n$ vektorene $\boldsymbol{\psi_m}$ er ortogonale og utgjør en basis for rommet med reelle (og komplekse) vektorer med $n$ elementer.
Parallelt med dette sjekker vi også det største (absolutte) avviket mellom denne matrisen og identitetsmatrisen, som er en annen enkel måte å gjøre sammenligningen på.
m = 1
w = 1
x, dx = np.linspace(-20, +20, 400, retstep=True)
v_infwell = 0*x
v_harmosc = 1/2*m*w**2*x**2
v_finwell = np.piecewise(x, [np.abs(x) > 2, np.abs(x) <= 2], [0, -10])
v_weird = np.random.rand(len(x))
potentials = [v_infwell, v_harmosc, v_finwell, v_weird]
titles = ["Infinite well", "Harmonic oscillator", "Finite well", "Weird potential"]
fig, axes = plt.subplots(1, 4, sharey=True)
fig.suptitle("Results of $\int \psi_m \psi_n \mathrm{d}x$")
axes[0].set_ylabel("n")
for i, v in enumerate(potentials):
_, waves = get_stationary_states(x, v, m)
mat = waves.dot(waves.T) * dx
# graphical presentation
cax = axes.flat[i].matshow(mat)
axes[i].set_title(titles[i])
axes[i].set_xlabel("m")
axes[i].xaxis.set_ticks([0, len(x)])
axes[i].yaxis.set_ticks([0, len(x)])
axes[i].xaxis.set_ticks_position('bottom')
# textual presentation
max_abs_dev = np.max(np.abs(mat - np.eye(len(x))))
print(f"Maximum deviation from identity matrix ({titles[i]}): {max_abs_dev}")
fig.colorbar(cax, ax=axes.ravel().tolist(), orientation="horizontal", pad=0.2);
Nå som vi har bekreftet ortogonaliteten til de stasjonære tilstandene, er vi godt forberedt til å beregne koeffisientene $c_n$ som behøves for å representere en vilkårlig starttilstand $\boldsymbol{\Psi_0}$ som en lineærkombinasjon $\boldsymbol{\Psi_0} = \sum_n c_n \boldsymbol{\psi_n}$ av de stasjonære tilstandene. For å finne $\boldsymbol{\Psi}$ for $t > 0$, trenger vi så bare sette på eksponentialfaktorene $e^{-i E_n t / \hbar}$.
Skriv en funksjon som beregner og returnerer koeffisientene $c_n$ som behøves for å representere en starttilstand $\boldsymbol{\Psi_0}$ som en lineærkombinasjon av de stasjonære tilstandene $\boldsymbol{\psi_n}$ i et diskretisert rom $\boldsymbol{x}$.
Skriv også en funksjon som beregner $\boldsymbol{\Psi}$ for en vilkårlig tid $t>0$ ut fra koeffisientene $c_n$, de stasjonære tilstandene $\boldsymbol{\psi_n}$ og energiene $E_n$. Merk at beregningene involverer tallet $i = \sqrt{-1}$, som i Python er representert med 1j.
Du kan godt verifisere at funksjonene i det minste klarer å reprodusere en eller annen starttilstand i et eller annet potensial ved $t = 0$. I de neste oppgavene vil du få testet dem også for $t > 0$.
Vi vet at $\boldsymbol{\psi_n}$ er ortogonale, dvs $\boldsymbol{\psi_m} \cdot \boldsymbol{\psi_n} = \sum_i \psi_m(x_i) \psi_n(x_i) \Delta x = \delta_{m n}$. Inspirert av innledningen "plukker vi ut" koeffisientene med $$\boldsymbol{\psi_m} \cdot \boldsymbol{\Psi_0} = \sum_n c_n \boldsymbol{\psi_m} \cdot \boldsymbol{\psi_n} = \sum_n c_n \delta_{m n} = c_m$$
For å beregne $\boldsymbol{\Psi}$ ved tiden $t$, beregner vi rett og slett bare $$\boldsymbol{\Psi} = \sum_n c_n \boldsymbol{\psi_n} e^{-i E_n t / \hbar}$$
Vi finner de stasjonære tilstandene og energiene i en harmonisk oscillator. Vi lager oss så en (rar) starttilstand $\boldsymbol{\Psi_0}$, beregner koeffisientene $c_n$ og bruker disse til å rekonstruere starttilstanden ved $t = 0$. Det gjenstår dermed bare å se at funksjonene fungerer også for $t > 0$ i de senere oppgavene.
def wave_coeffs(x, waves, wave):
dx = x[1] - x[0]
coeffs = waves.dot(wave) * dx
return coeffs
def wave_time(coeffs, waves, energies, t):
wave = waves.T.dot(coeffs * np.exp(-1j*energies*t/hbar))
return wave
m = 1
w = 1
x = np.linspace(-20, +20, 500)
v = 1/2*m*w**2*x**2
energies, waves = get_stationary_states(x, v, m)
wave0 = (3 - 4*x**2 + 3*x) * 5 * np.sin(0.5*x) * np.exp(1j * x)
coeffs = wave_coeffs(x, waves, wave0)
wave_reconstructed = wave_time(coeffs, waves, energies, 0)
fig, ax = plt.subplots()
ax.set_xlabel("$x$")
ax.set_ylabel("$\Psi(x)$")
ax.plot(x, np.real(wave0), linewidth=4, label="Original (real part)")
ax.plot(x, np.real(wave_reconstructed), linewidth=1, label="Reconstructed (real part)")
ax.plot(x, np.imag(wave0), linewidth=4, label="Original (imaginary part)")
ax.plot(x, np.imag(wave_reconstructed), linewidth=1, label="Reconstructed (imaginary part)")
ax.legend();
Nå har vi alt vi trenger for å studere tidsutviklingen av bølgefunksjonen.
Dette skal vi gjøre ved å lage en animasjon.
matplotlib er ikke et spesielt velegnet bibliotek til dette, men fungerer allikevel godt nok for vårt formål.
Nedenfor har vi satt sammen funksjonen animate_wave() som du skal modifisere til å gjøre nettopp dette.
Den tar inn et diskretisert rom $\boldsymbol{x}$, potensialverdier $\boldsymbol{V}$, en masse $m$, en starttilstand $\boldsymbol{\Psi_0}$ og antall bilder som skal vises per sekund, fps.
Funksjonen kaller animate() én gang per bilde i animasjonen, der en (foreløpig meningsløs) "bølgefunksjon" beregnes og vises i animasjonen.
Nedenfor har vi satt sammen litt mer informasjon om funksjonen og noen tips om hvordan du kan bruke den til å få tilfredsstillende resultater.
- Du kan vise real- eller imaginærdelen til funksjonen som animeres med argumentene
realogimag. - Du kan bestemme når animasjonen skal starte og slutte med argumentene
t1ogt2. Om du ikke ber om annet, vil den begynne ved $t_1 = 0$ og fortsette for evig og alltid til $t_2 = \infty$. Differansen mellom disse tallene bestemmer også varigheten på animasjonen i sekunder. - Du kan bestemme om animasjonen skal vises i sanntid eller forhåndsberegnes med argumentet
realtime.- I sanntidsmodus gjøres beregningene parallelt med animeringen. Du vil se resultater fra første sekund, men animasjonen kan bli hakkete hvis tidssteget er så lite at beregningene ikke klarer å holde følge med animeringen.
- I forhåndsberegningsmodus bruker maskinen tiden den trenger til å sette sammen alle bildene i animasjonen til en video. Tidssteget vil respekteres, men forhåndsberegningen kan ta en stund og du må bruke et endelig tidsintervall. Kjør animasjonene i forhåndsberegningsmodus før du leverer oppgaven, slik at de kan spilles av som en video.
Det er også noen tekniske ting som er greit å merke seg, i tilfelle ting ikke fungerer.
- Funksjonen
animate_wave()må kalles nederst i en Jupyter Notebook-celle for at animasjonen skal vises. - Hvis du irriterer deg over at figuren som animeres vises to ganger, kan du kommentere ut kodeblokken med
clear_output(). Vær i så fall klar over at dette vil fjerne all output som er generert av en celle fram tilclear_output()kalles. - Sanntidsmodus
realtime=Truekan kun brukes om du har bedt Jupyter Notebook eller Jupyter Lab bruke en interaktiv backend med%matplotlib notebookeller%matplotlib widgetøverst i notebooken. Om ikke må du forhåndsberegne animasjonene medrealtime=False.
Sett deg inn i og modifiser funksjonen animate_wave(), slik at den animerer $\boldsymbol{\Psi}$ fra starttilstanden $\boldsymbol{\Psi_0}$ i potensialet $V$ i rommet $\boldsymbol{x}$.
Vi beregner de stasjonære tilstandene $\boldsymbol{\psi_n}$, energiene $E_n$ og koeffisientene $c_n$ som representerer starttilstanden før animeringen begynner, da ingen av disse verdiene endrer seg under tidsutviklingen. Dette bruker vi så til å beregne $\boldsymbol{\Psi}$ ved hvert animasjonssteg. Vi plotter potensialet $\boldsymbol{V}$ og tilstanden $\boldsymbol{\Psi}$ i samme figur.
def animate_wave(x, v, m, wave0, fps, t1=0, t2=None, realtime=True, real=False, imag=False):
assert not (not realtime and t2 is None), "non-realtime animation must be finite in time"
dt = 1 / fps
nframes = None if t2 is None else int((t2 - t1) / dt) # None animates forever
nframesstr = "infinite" if t2 is None else f"{nframes}"
durationstr = "infinite" if t2 is None else f"{t2-t1:.2f}"
print("Animation information:")
print(f" Frames : {nframesstr}")
print(f" Framerate: {fps} FPS")
print(f" Duration : {durationstr}")
print(f" Time step: {dt}")
fig, ax = plt.subplots()
# prevent duplicate animation figure in non-realtime mode
#if not realtime:
#clear_output()
ax.set_xlabel("$x$")
ax.set_ylabel("$|\Psi|$, $\Re{(\Psi)}$, $\Im{(\Psi)}$")
# the commas matter!
ymax = max(np.max(np.abs(wave0)**2), np.max(np.abs(wave0)))
graph, = ax.plot([x[0], x[-1]], [0, +2*ymax]) # plot 2x wave0 to make room
if real:
graph2, = ax.plot([x[0], x[-1]], [0, -2*ymax]) # make room for negative values
if imag:
graph3, = ax.plot([x[0], x[-1]], [0, -2*ymax]) # make room for negative values
# plot extended potential with indication of infinite walls at ends
ax2 = ax.twinx()
v_max = np.min(v) + 1.1 * (np.max(v) - np.min(v)) + 1 # + 1 if v = const
x_ext = np.concatenate(([x[0]], x, [x[-1]]))
v_ext = np.concatenate(([v_max], v, [v_max]))
ax2.set_ylabel("$V(x)$")
ax2.plot(x_ext, v_ext, linewidth=3, color="black", label="V")
ax2.legend(loc="upper right")
energies, waves = get_stationary_states(x, v, m)
coeffs = wave_coeffs(x, waves, wave0)
def animate(i):
time = t1 + i*dt
wave = wave_time(coeffs, waves, energies, time)
# set_data() etc. modifies existing an existing object in a figure
# it is much more efficient than creating a new figure for every animation frame
graph.set_data(x, np.abs(wave))
graph.set_label(f"$|\Psi(x, t = {time:.2f})|$")
if real:
graph2.set_data(x, np.real(wave))
graph2.set_label(f"$\Re(\Psi(x, t = {time:.2f}))$")
if imag:
graph3.set_data(x, np.imag(wave))
graph3.set_label(f"$\Im(\Psi(x, t = {time:.2f}))$")
ax.legend(loc="upper left")
ani = matplotlib.animation.FuncAnimation(fig, animate, frames=nframes, interval=dt*1000, repeat=False)
if realtime:
return ani
else:
return HTML(ani.to_jshtml())
I resten av øvingen skal vi rett og slett bare teste ut animasjonsfunksjonen på en del forskjellige situasjoner.
Partikkel i boks¶
I regneøvingene har vi sett på superposisjonen av grunntilstanden og første eksiterte tilstand $$ \Psi(x, 0) = \frac{1}{\sqrt{2}} (\psi_1(x) + \psi_2(x)) $$ for en partikkel i boks. Vi fant så at sannsynlighetstettheten $|\Psi(x, t)|^2$ skulle oscillere fram og tilbake med perioden $$ T = \frac{2 \pi \hbar}{E_2 - E_1} $$
Preparér denne starttilstanden for et elektron i en boks og bekreft at perioden er riktig. Du kan for eksempel variere størrelsen på boksen for å finne en periode du har tid til å observere.
Vi lager en boks med en bredde som er liten nok til at (tallverdien til) perioden bare er på noen få sekunder. Så henter vi ut de stasjonære tilstandene og energiene til boksen, lineærkombinerer dem som angitt og setter i gang tidsutvikleren. Vi regner også ut den forventede perioden ved hjelp av energiene til tilstandene. Det ser ut til at bølgefunksjonen gjenvinner sin opprinnelige form etter den beregnede perioden!
x = np.linspace(0, 4, 400)
v = 0*x
m = 1
energies, waves = get_stationary_states(x, v, m)
wave = 1/np.sqrt(2) * (waves[0] + waves[1])
period = 2*np.pi*hbar/(energies[1] - energies[0])
print(f"Period: {period}")
animate_wave(x, v, m, wave, 10, t2=2*period, realtime=False)
Fri partikkel representert ved Gaussisk bølgepakke¶
I regneøvingene har vi også sett hvordan vi kan representere en fri partikkel som en gaussisk bølgepakke $$ \Psi(x, 0) = (2 \pi \Delta x^2)^{-1/4} e^{-(x-x_0)^2/4 \Delta x^2} e^{i p_0 x / \hbar} $$ med $$ \langle x \rangle = x_0 \quad \text{og} \quad \langle p \rangle = p_0 $$ Slik oppnådde vi også et "beste kompromiss" mellom usikkerhetene i posisjon $\Delta x$ og impuls $\Delta p$, nemlig $$ \Delta x \Delta p = \hbar / 2 $$ Dette er altså det nærmeste vi kommer å kunne representere en fri partikkel med en noenlunde veldefinert posisjon og impuls.
Plassér et elektron i en uendelig brønn med en gaussisk bølgepakke. Sørg for at den har plass til å bevege seg et stykke før den kræsjer i veggen. Studér utviklingen over tid. Hva skjer med formen til bølgepakken over tid? Hva skjer ved veggene? Hvilken hastighet har bølgepakkens tyngdepunkt?
Bølgepakken brer seg ut over tid. Dette svarer til fenomenet dispersjon. Bølgepakken er en kombinasjon av flere planbølger, hver med skarp impuls, så grovt sagt beveger ulike komponenter av totalbølgen seg med forskjellig hastighet. Dermed "diffunderer" bølgen ut over tid.
Denne tankegangen kan også overføres fysisk. I starten vet vi med stor nøyaktighet hvor partikkelen er, men vi vet mindre om impulsen og "hastigheten". Dermed blir vi mindre sikre på hvor partikkelen befinner seg over tid.
Ved veggene ser det ut til at det uendelige potensialet gjør at partikkelen spretter tilbake i motsatt retning som i en elastisk kollisjon.
Fra animasjonen ser vi at bølgepakken ser ut til å bevege seg med gruppehastigheten $\langle p \rangle / m$, som vi forventer fra klassiske betraktninger.
def gaussian(x, x0, stdx, p0):
return (2*np.pi*stdx**2)**(-1/4) * np.exp(-(x-x0)**2/(4*stdx**2)) * np.exp(1j * p0*x/hbar)
m = 1
x = np.linspace(-5, +45, 400)
v = 0*x
p0 = 3
wave = gaussian(x, 0, 1, p0)
animate_wave(x, v, m, wave, 10, t2=1.5*(x[-1]-x[0])/(p0/m), real=True, realtime=False)
