QM eines Teilchens in einer Dimension

David Gross

2021-05-04

In [1]:

import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg') 
plt.rcParams['figure.figsize'] = (9,6)


# Helpful for animations
# https://stackoverflow.com/questions/35532498/animation-in-ipython-notebook/46878531#46878531

In [27]:

# disable automatic display of plots
%matplotlib

import matplotlib.pyplot as plt
from matplotlib import animation, rc
import numpy as np
from scipy import linalg
from IPython.display import HTML, display, clear_output

def animateqm(psi,H,T,N,online=True):
    plt.rcParams['figure.figsize'] = (6,4)
    plt.rcParams['figure.dpi']=100
    plt.rcParams['animation.bitrate']=4000

    fig, ax = plt.subplots()
    ax.set_ylim((0,1))
    x = np.arange(0, len(psi))
    line, = ax.plot(x, np.abs(psi)**2)
    ts = np.linspace(0,T,N)

    def animate(t):
        evolved = np.abs(linalg.expm(-1j*t*H) @ psi)**2
        line.set_ydata(evolved)

    # Show individual frames while rendering.
    # Good for online demonstration, not suitable for web export
    if online:
        for t in ts:
            animate(t)
            display(fig)
            clear_output(wait = True)
    else: 
        anim = animation.FuncAnimation(fig, animate,
                                   frames=ts, interval=50, 
                                   blit=False) 
        return(HTML(anim.to_html5_video()))

Using matplotlib backend: agg

In [28]:

def laplace(n):
    return(
        np.array(
            [[ 1 if np.abs(np.mod(i-j,n)) in (1,n-1) else (-2 if i==j else 0) for i in range(n) ] for j in range(n) ]
        )
    )

# Center mu, variance s, momentum k
def wavePacket(mu,s,k,n):
    return(np.array([np.exp(-(x-mu)**2 / (2*s) + (2*np.pi*1j/n) * k * x) for x in range(n)]))

In [29]:

# free particle

n = 200;

V = np.array([0 for x in range(n)])

H = laplace(n) + np.diag(V);

psi = wavePacket(n/4,n/2,-10,n)

display(animateqm(psi, H, 300, 100, False))

In [30]:

# scattering at potential barrier 

n=500;

V=np.array([0 if x <= n/2  else .05 for x in range(n)])

H=laplace(n) + np.diag(V);

psi=wavePacket(n/2-60,30,-20,n);

display(animateqm(psi,H,300,100,False))

In [188]:

# tunneling

n = 500;

V = np.array([0 if (x <= n/2 or x >= n/2+2)  else 1 for x in range(n)])

H = laplace(n) + np.diag(V);

psi = wavePacket(n/2-60,30,-20,n);

display(animateqm(psi, H, 250, 100, False))