Quantenmechanik/ Numerik-Beispiel

Zurück zu Naeherungen | Hoch zum Quantenmechanik Anker | Vor zu Potenzial-Streuung

Dieser Abschnitt tastet sich numerisch an die Schrödinger-Gleichung heran in einer Dimension und findet Lösungen im Potenzialtopf. Das Beispiel Morse-Potenzial lässt sich zum besseren Vergleich auch analytisch lösen. Ein Python-Skript soll zeigen, dass mit rustikaler Numerik recht gute Ergebnisse herauskommen -- bei niedrigen Energien mit der eindimensionalen Schrödinger-Gleichung, wo z.B. die WKB-Näherung zu schlecht wäre.

Ein radiales Potenzialmodell für die Anziehung zwischen Atomen oder Molekülen mit einem abstoßendem Kernbereich ist das Lennard-Jones-Potenzial mit den zwei Parametern E,s:

${\displaystyle V_{LJ}(r)\;:\quad V_{LJ}(r)/E=4[(s/r)^{12}-(s/r)^{6}]}$

Es ist negativ für r>s und positiv für r<s.
Sein Minimum (-E) liegt bei ${\displaystyle r=2^{1/6}s.}$

Hier soll mit der ähnlichen Morse-Potenzialfunktion gerechnet werden:

${\displaystyle p(x)=\exp(-2x/b)-2\exp(-x/b)}$.

Schrödinger-Gleichung: ${\displaystyle (-\hbar ^{2}/(2m))\cdot \phi ''+Vp(x)\phi =E\phi }$.

Der Potenzialtopf ist sehr unsymmetrisch, er steigt bei negativem ${\displaystyle x<(-b)}$ steil an, hat ein Minimum (-V) bei x=0 und flacht für großes ${\displaystyle x>b}$ gegen Null ab. Die Länge b misst die Ausdehnung des Topfes, die Konstante V seine Tiefe. Ein solches Potenzial simuliert in etwa die Kraft, mit der zwei Atomkerne in einem Molekül gebunden sind. Es kann ein typisches Spektrum der Vibrationen des Moleküls darstellen.

Die X-Variable wird zu dimensionslosem y skaliert:

${\displaystyle y=x/b;\;d/dy=b(d/dx);\;p(y)=\exp(-2y)-2\exp(-y)}$.

Gleichung: ${\displaystyle [-\hbar ^{2}b^{2}/(2m)]\phi ''(y)+Vp(y)\phi =E\phi }$.

Energie-Skalierung, ebenfalls zu dimensionslosen Werten:

${\displaystyle ({\bar {E}},{\bar {V}})=(E,V)\cdot 2m/(\hbar ^{2}b^{2})}$

Ergibt einfach, mit unterdrückten Strichen: ${\displaystyle -\phi ''=(E-Vp(y))\phi }$.

Das Potenzial-Minimum ist (-V) bei y=0. Laut Theorie sind die Niveaus:

${\displaystyle E_{n}=-({\sqrt {V}}-(n+{\tfrac {1}{2}}))^{2};\quad {\sqrt {-E_{n}}}=[{\sqrt {V}}-(n+{\tfrac {1}{2}})]\quad (n=0...N)}$

Mit dem Beispiel V=200 ergeben sich 14 Energie-Niveaus n=0...13.

Das Skript löst die Schrödinger-Gleichung im Morse-Potenzial und findet eine Serie der ersten 8 gebundenen Zustände. Die Wellenfunktionen und die Energieniveaus werden als SVG-Grafik ausgegeben. Trotz der langsamen Sprache Python und der lahmen Technik von numerischer Integration tausender Probeschüsse ist das Ergebnis in ein paar Sekunden fertig. Im ersten Durchgang wird für ein Raster von 50 Punkten der Energie im Potenzialtopf je eine Probewelle berechnet.

Angenäherte Differenzgleichung: die Gleichung ${\displaystyle f''=g(x,f)}$ wird numerisch auf einem Gitter von gleich entfernten x-Punkten zur Differenzgleichung

${\displaystyle (f[x+d]+f[x-d]-2f[x])/d^{2}=g(x,f[x])}$.

Sind zwei Anfangswerte bekannt, kann eine Kurve f() Schritt für Schritt 'integriert' werden:

${\displaystyle f[x+d]=d^{2}\cdot g(x,f[x])+(2f[x]-f[x-d])}$.

Probewelle: An den Rändern des Topfes, das heißt an den klassischen Umkehrpunkten, wird zuerst je eine (mindestens exponentiell) abfallende Lösung in die Potenzialmauer hinein iteriert. Dann wird die oszillierende Welle f von Rand zu Rand geführt, mit stetig-differenzierbarer Anfangsbedingung links. Am rechten Punkt wird der Sprung der relativen Steigung ${\displaystyle f'(x)/f(x)}$ bezüglich der Randlösung vermerkt. Wäre die Welle ein erlaubter Zustand, dann wäre dieser Fehler Null. Gleichzeitig werden auch die Knoten (Nulldurchgänge) der Welle im Topf gezählt.

Im zweiten Durchlauf werden benachbarte Paare aus dem Raster gezogen, die zur gleichen Knotenzahl gehüren, aber zwischen denen sich das Vorzeichen beim Steigungs-Knacks umkehrt. Ein solches Paar von Energie-Niveaus ist dann eine Klammer, zwischen der sich ein erlaubtes Niveau mit Knotenzahl N befindet. Mit Dichotomie und einer Iteration von Probewellen wird die Eigenfunktion und der Eigenwert der Energie ermittelt. Genauer, die Nullstelle der Steigungs- Unstetigkeit. Spontan ergeben sich N=0,1...7 als Serie von Quantenzuständen.

Die Fehler der Energie, verglichen mit der Theorie, sind kleiner als 0,1%. Die höchsten Niveaus N=8 bis 13 werden nicht erkannt. Das Raster nahe beim Niveau Null müsste dazu feiner gewählt werden. Besser noch wäre, beide Durchläufe zu vereinen und das Such-Raster adaptiv nach den Energie-Abständen zu justieren. Vielleicht adaptive Schrittweiten auch bei den unzähligen Summierungen der Differenzengleichung, das Näherungsverfahren für die Integration der Differenzialgleichung?

Löse die Gleichung: ${\displaystyle -\phi ''(z)+Vu(z)\phi =E\phi ;\;u(z)=e^{-2z}-2e^{-z}}$.

Die Strategie: schalte um zu bequemeren Variablen, im Definitionsbereich von z zu t und im Wertebereich von ${\displaystyle \phi (z)}$ zu ${\displaystyle w(t)}$. Die transformierte Differenzialgleichung in {w,t} sollte lösbar sein.

Ansatz: ${\displaystyle t=t(z)=fe^{-z}\quad {\text{und}}\quad w(t)=e^{at}t^{b}\phi (z)}$.

Die drei freien Parameter ${\displaystyle a,b,f}$ werden später getrimmt. Ableitungen nach z werden mit Strichen, solche nach t mit Punkten markiert.

${\displaystyle (dt/dz)=-fe^{-z}=-t;\;(d^{2}t/dz^{2})=fe^{-z}=t;\;e^{-z}=(t/f);\;z=-\ln(t/f)}$

${\displaystyle w'=(dw/dz)=at'w+bw(t'/t)+w(\phi '/\phi )=(-at-b)w+(\phi '/\phi )w}$

${\displaystyle w''=[-at-b+(\phi '/\phi )]w'+[-at'+(\phi ''\phi -\phi '^{2})/\phi ^{2}]w}$

(vorgemerkt)${\displaystyle \quad w'=-{\dot {w}}t;\;w''=(d/dz)(-{\dot {w}}t)={\ddot {w}}t^{2}+{\dot {w}}t}$

Mit ${\displaystyle (\phi '/\phi )=(w'/w)+(at+b)\quad {\text{und}}\quad (\phi ''/\phi )=Vu(t)-E=:P(t)}$, das ein Polynom der Ordnung 2 in t ist:

${\displaystyle w''=w'^{2}/w+w[at+P(t)-((w'/w)+(at+b))^{2}]}$

${\displaystyle w''=({\dot {w}}^{2}t^{2})/w+w[at+P(t)-((-{\dot {w}}t/w)+(at+b))^{2}]}$

${\displaystyle w''=w[at+P(t)-(at+b)^{2}+2(at+b)t({\dot {w}}/w)],\quad }$ Term in ${\displaystyle {\dot {w}}^{2}}$ entfällt.

Mit ${\displaystyle w''={\ddot {w}}t^{2}+{\dot {w}}t}$ wird schließlich die Gleichung sauber in w und t:

${\displaystyle {\ddot {w}}t^{2}+{\dot {w}}[t-2(at+b)t]+w[(at+b)^{2}-at-P(t)]=0}$

Nun ist ${\displaystyle P(t)=V((t/f)^{2}-2(t/f))-E}$. Mit der Wahl ${\displaystyle b^{2}+E=0}$ kann der inhomogene Term entfallen und der Rest durch t geteilt werden. Daher bleibt:

${\displaystyle {\ddot {w}}t+{\dot {w}}[1-2(at+b)]+w[(at+v)^{2}-at-Vt/f^{2}+2V/f]=0}$

Es gibt Licht am Ende des Tunnels.

Weitere Wahl: ${\displaystyle a^{2}-V/f^{2}=0}$ entfernt den Faktor von ${\displaystyle t^{2}}$.

Dritte Wahl : ${\displaystyle 2a=1}$ entfernt auch den Faktor von t im w-Term. Also a=1/2.

${\displaystyle {\ddot {w}}t+{\dot {w}}[1+b-t]+w[b-a+2fa^{2}]=0;\;f=\pm 2{\sqrt {V}}}$.

Was hier bleibt, ist die Differenzialgleichung der Kummer-Funktionen F.

${\displaystyle F(A,B;z)\;}$ löst die Gleichung ${\displaystyle \;zF''+(B-z)F'-AF=0}$.

Wenn dazu noch A eine nicht-positive Ganzzahl (-n) ist, ist die Lösung F ein Polynom in z.

Im Abschnitt Potenzial-Streuung steht etwas mehr zu den Kummer-Funktionen.

In der Hoffnung, dass mit den Polynomen zu quanteln ist, rechnet man (n) aus:

${\displaystyle f=\pm {\sqrt {V}}/a;\quad n=b-a-2fa^{2}=(\pm {\sqrt {-E}})-(1/2)+(\pm {\sqrt {V}})}$

Die Physik sagt: E ist eine negative Bindungs-Energie zwischen -V und 0.

Die sinnvollen Vorzeichen sind daher: ${\displaystyle b=-{\sqrt {-E}}\quad {\text{und}}\quad f=2{\sqrt {V}}}$.

Erlaubte Energie-Niveaus: ${\displaystyle {\sqrt {-E_{n}}}={\sqrt {V}}-(n+{\tfrac {1}{2}});\;n=0,1,2...}$

Die Energien sind negativ und die Serie der ${\displaystyle E_{n}}$ ist endlich.

Sind die Wellenfunktionen ${\displaystyle \phi (z)=w(t)e^{-at}t^{-b}}$ nun überhaupt gutartig, wobei ${\displaystyle w(t)=F(t)}$ ein Polynom ist?

${\displaystyle t=fe^{-z};\;z=-log(t/f);\;a=1/2;\;b=-{\sqrt {-E}};\;f=2{\sqrt {V}}>0}$.

${\displaystyle \phi (z)=F(t)e^{-t/2}t^{\sqrt {-E}}\;}$ mit positiver Potenz von t.

Für ${\displaystyle z\to \infty ,\;t\to 0\;\Rightarrow \phi \to F(t)t^{\sqrt {-E}}=0}$

Für ${\displaystyle z\to -\infty ,\;t\to \infty \;\Rightarrow \phi \propto e^{-t/2}\to 0}$.

Tatsächlich sind annehmbare zeitunabhängige Schrödinger-Lösungen gefunden.

from math import exp,sqrt,log
def div(p,q) : return int(p/q)
def pot(x) : return exp(-2*x)-2*exp(-x)

def shoot(x,dx,f,df,xlim, v0,erg) : # decay tails in both dirs
  xa=x;xb=x+dx;xc=xb+dx; ya=f;yb=f+df; steps=0; smax=1000; ymax=3*f; node=0
  smax= int(abs((xlim-x)/dx)); ix=0; done=False; zb=abs(yb); data=[xa,ya]
  while not done : # f"= (e-u)f = (yc+ya-2yb)/dx^2 = (erg-u(xb))*yb
    z= -dx*dx*yb*(erg-v0*pot(xb)); yc= z+2*yb-ya;xc=xb+dx; zc=abs(yc); steps+=1
    if (yc*yb)<0 : node+=1
    done= (steps>smax)or(zc>zb) or (yc<0)or(yc>ymax); ix+=1
    if not done : xa=xb;ya=yb; xb=xc;yb=yc; zb=zc; data += [xa,ya] 
    if done and (steps<smax) : yc= yc*smax/(steps+0.0)  
  return yc,node,ix, data 

def waveshot(x,f,df, y, v0,erg) : # start f(x), df=derivative, oscillates
  steps=0; smax=1000; node=0; dx=(y-x)/(smax+0.0); dfdx=df*dx
  xa=x-dx;xb=x; ya=f-dfdx; yb=f; ix=0; noisy=False; data=[xb,yb]
  while (ix<smax) : # f"= (u-e)f = (yc+ya-2yb)/dx^2 = (u(xb)-erg)*yb
    z= dx*dx*yb*(v0*pot(xb)-erg); yc= z+2*yb-ya; xc=xb+dx; ix +=1
    if (yc*yb)<0 : node+=1
    xa=xb;ya=yb; xb=xc;yb=yc; data += [xb,yb]
    if noisy and (ix%50)==0 : print('z='+str(z)) # should be wavy
  return xb,yb,(yb-ya),dx, node, data 

def bestshot(x,dx,f, v0,erg) : #  try to get optimal decay to zero
  df=0; targ=0; dfhi=0; dflo=-f; noisy=0 # zhi=2*f; zlo=-0.1*f; 
  zbefore=10; i=0; busy=True; zmi=targ-0.1; zmx=targ+0.1;zhi=zmx;zlo=zmi
  while busy :
    z,node,ix,data = shoot(x,dx,f,df, 5.0, v0,erg);
    if noisy>=2 : print(str([z,node,ix,df]))
    if z>targ : dfhi=df;zhi=z 
    else : dflo=df;zlo=z
    df=0.5*(dfhi+dflo);
    if (zlo>zmi)and(zlo<zhi)and(zhi<zmx) :
      df=dflo+(targ-zlo)*(dfhi-dflo)/(zhi-zlo)
    i+=1; busy= (abs(z-zbefore)>1e-8); zbefore=z
  if noisy>0 : print('Shotloop i,df,dx,z='+str([i,df,dx,z])) 
  return df,dx, data

def schroedinger(erg,v0): # data 3 arrays of xy pairs, queue-wave-queue.
  dx=0.01; xa=0; xb=0; fa=1; fb=1; noisy=False; data=[[]]*3
  while erg> (v0*pot(xa)) : xa += dx
  while erg> (v0*pot(xb)) : xb -= dx 
  if noisy : print('erg,v0,xb,xa='+str([erg,v0,xb,xa]))
  dfa,dxa,data[2] = bestshot(xa,dx,fa, v0,erg);
  dfb,dxb,data[0] = bestshot(xb,-dx,fb, v0,erg); 
  xx,yy,dy,dxx, node, data[1] = waveshot(xb,fb,dfb/dxb, xa, v0,erg)
  if noisy : print('waveshot to '+str([xx,yy,dy,node]))
  ta=dfa/dxa/fa; ht=dy/dxx/yy # compare target to hit
  # print('Target and hit slopes: '+str([ta,ht]))
  return node, ht/ta, data

def pairflip(d) :
  j=len(d); k=j-2;i=0
  while (i<k) :
    x=d[i];y=d[i+1];d[i]=d[k];d[i+1]=d[k+1];d[k]=x;d[k+1]=y
    i=i+2; k=k-2 

def squareint(d) : # square-integral
  j=len(d); i=0; sum=0
  while i<(j-3) :
    dx= d[i+2]-d[i]; y=d[i+1];z=d[i+3]; sum += 0.5*dx*(y*y+z*z); i+=2
  return sum

def rescale(d,f) : # return bounds too
  j=len(d); i=0; xmin=1e12;xmax=-xmin;ymin=xmin;ymax=xmax
  while i<j :
    x=d[i]; y= d[i+1]*f; d[i+1]=y; i=i+2
    xmin= x if(x<xmin) else xmin; xmax= x if (x>xmax) else xmax
    ymin= y if(y<ymin) else ymin; ymax= y if (y>ymax) else ymax
  return xmin,ymin, xmax,ymax

def findwave(v,scan,ifd,node,plot) : #bracket search for solution in scan table
  ea,na,ra= tuple(scan[ifd]); eb,nb,rb= tuple(scan[ifd+1])
  ra= ra-1; rb=rb-1; k=0; ec=0.5*(ea+eb); r=1; ed=1
  while (k<20)and(abs(r)>1e-5) :
    print('k, r,ea,eb,ec,ed='+str([k,r,ea,eb,ec,ed]))
    nn,r,data = schroedinger(ec,v); r= r-1; k+=1 # r seems negative always!?
    if (r*ra)>0 : ra=r;ea=ec
    else : rb=r; eb=ec
    ed= ea+(-ra)*(eb-ea)/(rb-ra) # better target zero?
    if ((ed-ea)*(ed-eb))<0 : ec=ed
    else : ec= 0.5*(ea+eb)
  print('Data ranges:') # reverse first list, y-rescale last
  pairflip(data[0]); d=data[1]; j=len(d); yb=d[j-1]; sum=0
  rescale(data[2],yb)
  for k in range(len(data)) :
    d=data[k]; j= len(d); xa,ya,xb,yb=d[0],d[1],d[j-2],d[j-1]
    print('j xy xy'+str([j,xa,ya,xb,yb]))
    s=squareint(d); print('  SquareInt='+str(s)); sum+=s
  norm=sqrt(sum); print('Norm='+str(norm)) # divide all stuff by this one
  for k in range(3) :
    xa,ya,xb,yb= rescale(data[k],1.0/norm); print(str([xa,ya,xb,yb]))  
    if k==0 : xmin,ymin,xmax,ymax= xa,ya,xb,yb
    else :
      xmin=min(xmin,xa);ymin=min(ymin,ya);xmax=max(xmax,xb);ymax=max(ymax,yb)
  if plot : waveplot('morsew'+str(node), node, ec, xmin,ymin,xmax,ymax,data)
  return([ec, node, xmin,ymin,xmax,ymax,data])

def multiplot(d) :
  n=len(d); xmin=1e12;xmax=-xmin;ymin=xmin;ymax=xmax
  txt=[]; fn='morse'; ix=[0]*(1+3*n)
  for i in range(n) :
    j=1+3*i;k=i+1; ix[j],ix[j+1],ix[j+2]= k,k,k # color index
    ec, node, xa,ya,xb,yb,dt = tuple(d[i])
    xmin= xa if(xa<xmin) else xmin; xmax= xb if (xb>xmax) else xmax
    ymin= ya if(ya<ymin) else ymin; ymax= yb if (yb>ymax) else ymax
  xyf=[ xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin]; xys=[xyf]
  for i in range(n) :
    ec, node, xa,ya,xb,yb,dt = tuple(d[i]); xys+=[dt[0],dt[1],dt[2]]
    txt+= ['E'+str(node)+'= '+str(aprx(ec))]
  xs,fx,xpix= xmin, 750/(xmax-xmin),25 # scale user units, min in pixels
  ys,fy,ypix= ymin,-550/(ymax-ymin),575 ; scale=(xs,fx,xpix, ys,fy,ypix)
  colo=['#aaa','#00a','#077','#0a0','#770','#a00','#707','#0ff','#f0f']
  sd=svgdump(); buf=sd.plotbegin(800,600,1); i=0; dy=(ymax-ymin)/20.0
  for j in range(len(txt)) :
    buf += label(sd,scale, txt[j], 0.2*xmin+0.8*xmax,ymax-dy-j*dy)
  for t in xys : buf += curve(sd,scale, t, fill=colo[ix[i]]); i+=1
  g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()

def energytheory(v,d) : # compare energy, theory versus numerical
  n=0; ok=True; ec=0; ld=len(d)
  while ok : 
    if n<ld : ec, node, xa,ya,xb,yb,dt = tuple(d[n])
    z= sqrt(v)-(n+0.5); ok= z>=0
    if ok : y=-z*z; print(str([n,y,ec])); n+=1; ec=0

def test() : # stationary schroedinger equ numerics, 8 energy levels found
  v=200; n=50; scan=[[]]*n; waves=[] # scan 50 energy levels
  for i in range(1,n) : # scan for node numbers and slope defects
    erg= -v + i*v/(n+0.0); node,rsl, data= schroedinger(erg,v)
    scan[i]=[erg,node,rsl]
    print('*** Erg,nodes,rslope='+str(scan[i]))
  for i in range(1,n-1) :
    ea,na,ra= tuple(scan[i]); eb,nb,rb= tuple(scan[i+1])
    found= (na==nb)and ((ra-1.0)*(rb-1.0)<0) # cross sloperatio=1
    if found :
      print('found i='+str(i)+',nodes='+str(na)); ifd=i;node=na
      wave= findwave(v,scan,ifd,node, False); waves += [wave]
  print('Number of waves='+str(len(waves))+'. Theory vs. Numerical:');
  multiplot(waves); energytheory(v,waves)

###### plot

class svgdump :
  def plotbegin(sd, w,h,f) : # args: width height zoomfactor
    W=str(w); H=str(h); u=''; fw=str(f*w); fh=str(f*h)
    s=( '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n'+
     '<!DOCTYPE svg>\n'+
     '<svg xmlns="http://www.w3.org/2000/svg"\n'+
     '  xmlns:xlink="http://www.w3.org/1999/xlink"\n'+
     '  width="'+W+'" height="'+H+'" viewBox="0 0 '+fw+' '+fh+'">\n')
    return s

  def plotend(sd) : return '</svg>\n'
  def move(sd,x,y) : return 'M'+str(x)+','+str(y)
  def line(sd,x,y) : return 'L'+str(x)+','+str(y)
      
  def labl(sd,s,x,y, size=14,fill='#000') :
    f=' font-family="Arial"'
    t='<tspan x="'+str(x)+'" y="'+str(y)+'">'+s+'</tspan>'; fsz= str(size)
    return ('<text fill="'+fill+'"'+f+' font-size="'+fsz+'px">'+t+'</text>\n')

  def path(sd,w,s,lc='#000') : # path w=with lc=linecolor
    t='<path fill="none" stroke="'+lc+'" stroke-width="'+str(w)+'" '
    t+= ' stroke-linecap="round" ' 
    return t + 'd="\n'+s+'" />\n'

# end class svgdump

def aprx(x) : # 3-digit approximation
  y= (-x) if(x<0) else x; z= int(1e3*y+1e-3); z= (-z) if (x<0) else z
  return z/1000.0

def curve(sd,scale,xy, fill='#000', width=1) :
  (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= div(len(xy),2) 
  for i in range(n) : 
    x= aprx(xmin+fx*(xy[2*i]-xs)); y= aprx(ymin+fy*(xy[2*i+1]-ys))
    s+= sd.move(x,y) if(i==0) else sd.line(x,y)
  return sd.path(width,s,lc=fill)

def label(sd,scale,txt,tx,ty) : # tx,ty plot coordinates not pixels
  (xs,fx,xmin, ys,fy,ymin) = scale 
  x= aprx(xmin+fx*(tx-xs)); y= aprx(ymin+fy*(ty-ys))
  return sd.labl(txt,x,y)

######## main program #####
test()

### optional 

def waveplot(fn, lvl, erg, xa,ya,xb,yb,data) :
  txt=['Wavefunct Level '+str(lvl)+' Energy '+str(aprx(erg)),]
  xyf=[ xa,ya, xb,ya, xb,yb, xa,yb, xa,ya ] # frame
  xys=[xyf,data[0],data[1],data[2]]
  xs,fx,xpix= xa, 750/(xb-xa),25 # scale user units, min in pixels
  ys,fy,ypix= ya,-550/(yb-ya),575 ; scale=(xs,fx,xpix, ys,fy,ypix)
  colo=['#777','#00a','#077','#0a0','#770','#a00']
  sd=svgdump(); buf=sd.plotbegin(800,600,1); i=0; dy=(yb-ya)/20.0
  for j in range(len(txt)) :
    buf += label(sd,scale, txt[j], (xa+xb)*0.5,yb-dy-j*dy)
  # for j in range(len(lbl)) :
  #  buf += label(sd,scale, lbl[j], log(fg[j]),yb)
  for t in xys : buf += curve(sd,scale, t, fill=colo[i%7]); i +=1
  g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()

Beim Wasserstoff-Atom bilden die radialen Schrödinger-Lösungen im zentralen Coulomb-Potenzial Folgen um Laguerre-Polynome herum. Sie sind 'gequantelt' mit zwei ganzzahligen Indizes, worin die Hauptquantenzahl n und die Drehimpuls-Zahl l vorkommen. Wegen der besondere Symmetrie des Potenzials sind die stationären Niveaus stark entartet -- zu gegebenem n haben l=0,1,...n-1 die selbe Bindungsenergie.

Die Alkali-Atome Li bis Cs sind dem Wasserstoff darin ähnlich, dass ein einzelnes Elektron-Niveau sich von den anderen abhebt und in einem effektiven Zentralpotenzial der restlichen Elektronenwolken angenähert beschrieben werden kann. Das Leuchtelektron hat einen deutlichen Energie-Abstand zu den anderen Niveaus, die nach dem Pauli-Prinzip vorher dran kamen und eine stabile, wenig reaktionsfreudige Edelgas-Konfiguration belegen.

Ein Modell dieses effektiven radialen Potenzials soll kurz besprochen werden. Zur Anziehung durch den Atomkern kommt eine Abstoßung durch die inneren Elektronen. Die radiale Differenzialgleichung ist immer noch geschlossen lösbar, mit Hilfe von Modifizierten Laguerre-Polynomen. Die Entartung von ${\displaystyle (n,l)}$ Paaren wird aber aufgehoben.

Der Atomrumpf soll mit zwei dimensionslosen Parametern auskommen, nennen wir sie U und C. Statt der Kernladungszahl Z soll der effektive Coulomb-Anteil, also der in ${\displaystyle (1/r)}$, eine viel kleinere Anziehung (Z-U) haben. Der abstoßende Anteil soll mit ${\displaystyle (1/r^{2})}$ stark abfallen und einen empirischen Faktor C bekommen. Es ist zu prüfen, ob ein Paar ${\displaystyle \{U,C\}}$ beispielsweise einen Bereich des Natrium-Spektrums gut annähert oder ob die Werte noch feiner variiert werden müssen, ob also ein Fit mit mehr als 2 Parametern benötigt wird.

Zuerst mag man die Differenzialgleichung dimensionslos umschreiben, was allgemein für numerische Berechnungen Vorteile bringt. Die Länge wird skaliert in Bezug auf den Bohr-Radius ${\displaystyle a={\frac {\hbar ^{2}}{2m}}{\frac {4\pi \epsilon _{0}}{e^{2}}}}$. Die Energie wird ausgedrückt relativ zur Wasserstoff-Ionisierungsenergie

${\displaystyle E_{0}={\frac {\hbar ^{2}}{2ma^{2}}}}$. Damit ist ${\displaystyle {\frac {\hbar ^{2}}{2m}}=E_{0}a^{2}.\quad E_{0}=}$ 13,6 eV; a=0,0529 nm.

Vernachlässigt wird, die Vielteilchen-Gleichung in Schwerpunktskoordinaten und Relativkoordinaten von Rumpf und Leuchtelektron zu separieren und statt der Elektron-Masse eine reduzierte Masse einzusetzen. Vernachlässigt wird auch der Elektron-Spin, der spinorwertige Winkelfunktionen ${\displaystyle Y_{j;lm}}$ an Stelle der skalaren ${\displaystyle Y_{lm}}$ erfordern würde. Der Spin brächte darüber hinaus eine Feinstruktur ein mit Spin-Bahn-Kopplungsterm im Hamilton-Operator.

Potenzialmodell: ${\displaystyle V={\frac {-(Z-U)q^{2}}{4\pi \epsilon _{0}r}}+{\frac {\hbar ^{2}}{2m}}{\frac {C}{r^{2}}}={\frac {-(Z-U)\hbar ^{2}}{amr}}+{\frac {\hbar ^{2}C}{2mr^{2}}}}$

${\displaystyle {\frac {V}{E_{0}}}=-(Z-U){\frac {2a}{r}}+C{\frac {a^{2}}{r^{2}}};\quad r=|{\vec {x}}|}$

Energie-Niveaus ${\displaystyle F:=(E/E_{0})}$ kommen aus der Eigenwert-Gleichung

${\displaystyle \left[-\nabla ^{2}a^{2}-{\frac {(Z-U)(2a)}{r}}-{\frac {Ca^{2}}{r^{2}}}-F\right]\Psi =0}$

Mit skaliertem Radius ${\displaystyle s=(r/a);\;{\vec {s}}={\vec {x}}/a}$ und ${\displaystyle (d/dr)=a(d/ds)}$:

${\displaystyle \left[-\nabla _{s}^{2}-{\frac {2(Z-U)}{s}}-{\frac {C}{s^{2}}}-F\right]\Psi =0}$

Ansatz: ${\displaystyle \Psi =Y_{lm}(\theta ,\phi )u(s)}$.

Radialgleichung: ${\displaystyle [{\frac {d^{2}}{ds^{2}}}+{\frac {2}{s}}{\frac {d}{ds}}-{\frac {l(l+1)}{s^{2}}}]u(s)+[{\frac {2(Z-U)}{s}}-{\frac {C}{s^{2}}}+F]u(s)=0}$

Nach einem Muster, das schon beim Wasserstoff-Atom geholfen hat, wird noch einmal die Radius-variable skaliert ${\displaystyle x=s/b}$. Und von der Funktion u werden ein Potenz-Faktor und ein Exponential-Faktor abgespalten:

${\displaystyle u(s)=R(x)x^{\lambda }\exp(-x/\mu )}$.

Als beste Option kommt folgendes heraus, um die Eigenwerte F<0 einzukreisen:

${\displaystyle x=2s{\sqrt {-F/E_{0}}};\;u=R(x)x^{\lambda }\exp(-x/2)}$

Damit ergibt sich folgende Differenzialgleichung für ${\displaystyle R(x)}$:

${\displaystyle \left[x{\frac {d^{2}}{dx^{2}}}+(2\lambda +2-x){\frac {d}{dx}}+{\frac {\eta }{x}}+\nu \right]R(x)=0;\quad \eta :=\lambda (\lambda +1)-l(l+1)-C;\quad \nu :={\sqrt {-E_{0}/F}}(Z-U)-\lambda -1.}$

Dies ist fast die definierende Gleichung für die Kummer-Funktionen, siehe deren Definition im Abschnitt über Streuung.

${\displaystyle K(A,B;z)\;}$ löst die Gleichung ${\displaystyle \;zK''+(B-z)K'-AK=0}$.

Ist ${\displaystyle (-A)}$ eine ganze Zahl n=0,1,2,3... dann ist K ein modifiziertes Laguerre- Polynom.

Um solche Lösungen anzupeilen, wird gesetzt: ${\displaystyle \eta =0}$ und ${\displaystyle \nu =n=0,1,2,3...}$

• Für ${\displaystyle \lambda }$ folgt: ${\displaystyle \lambda ={\sqrt {(l+{\tfrac {1}{2}})^{2}+C}}-{\tfrac {1}{2}}}$.

Mit der positiven Wurzel wie hier ist ${\displaystyle u(s)}$ bei Null regulär und normierbar wegen des abfallenden Exponentialanteils. ${\displaystyle R(x)=R_{n\lambda }(x)=K(-n,2\lambda +2;x)}$ sind Polynome. Aus ${\displaystyle {\sqrt {-E_{0}/F}}=(\lambda +1+n)/(Z-U)}$ folgen

• Die Energie-Eigenwerte ${\displaystyle F_{n\lambda }=-E_{0}\left({\frac {Z-U}{n+\lambda +1}}\right)^{2}}$.

Beim Wasserstoff-Atom gibt es in der Energie-Reihe die Hauptquantenzahl ${\displaystyle N=n+l+1}$; in diesem Modell mit Abschirmung ist die effektive Zahl im Nenner ${\displaystyle N'=n+\lambda +1}$ nicht mehr ganz und im Vergleich größer. Nimmt man die Zahl N als Referenz, dann wird die Entartung der Drehimpulse l aufgehoben, die höheren Drehimpuls-Zahlen senken die Energie ab.

Das Natrium-Atom hat den Grundzustand 3s und bekannte angeregte Orbitale 3p 4s 3d 4p 4f. Diese werden versuchsweise mit dem Schema des Modells und passenden Parametern ${\displaystyle l,n,\lambda ,Z-U,C}$ konfrontiert. (Die Spin-Orbit-Aufspaltung in ${\displaystyle 3p_{1/2}}$ und ${\displaystyle 3p_{3/2}}$, die für die doppelte D-Linie im Spektrum verantwortlich ist, wird hier unter den Teppich gekehrt.)

In der ersten Zeile wurde wohl C erraten und (Z-U) gefittet, in den folgenden eher umgekehrt. Mit konstanten Werten U,C scheint es unmöglich, die Liste zu reproduzieren. Ein Zweiparameter-Modell ist doch etwas primitiv. Es nutzt der Didaktik, wegen des lösbaren Spektrums.

Das Beispielskript (dblslit.py) simuliert den Doppelspaltversuch auf einem rechteckigen Gebiet. Die zwei Raumdimensionen haben ein Raster von 120 mal 90 Punkten und die Zeit wird in 240 diskreten Schritten hochgezählt. Die Skriptsprache Python beherrscht komplexe Arithmetik sehr gut und ist schnell genug, um den ganzen Vorgang in etwa einer Minute auf einem Raspberry Pi durchzurechnen.

Die Schrödinger-Gleichung wird durch eine Differenzengleichung angenähert. Integration über die Zeit erfolgt im halb-impliziten Verfahren nach Crank und Nicholson, das die Unitarität ziemlich gut sichert, die Erhaltung der Summe aller Betragsquadrate. Für die zwei Raumdimensionen wird der Hamilton-Operator aufgespalten in den x-abhängigen und den y-abhängigen Teil. Jeder Zeitschritt hat zwei Teilschritte, in denen der eine oder der andere Halb-Operator das Fortschreiten der Wellenfunktion bestimmt. Es kommen daher nur eindimensionale Schrödinger-Differenzen-Gleichungen vor, und zwar als einfach zu lösende tridiagonale lineare Gleichungssysteme.

Das Spielfeld x=0...120, y=0...90 hat eine Potenzialmauer bei x=80, mit einem dreieckigen Profil, das 5 Punkte dick ist in der x-Richtung. Die Mauer hat zwei Löcher, durch die ein Teil der Welle austreten soll. Im Gebiet x<80 wird ein Wellenpaket mit Gaußkurven-Amplitude in zwei Dimensionen aufgesetzt, multipliziert mit dem Phasenfaktor exp(ikx). Der Wellenvektor k entpricht einer Wellenlänge von 20 Rasterpunkten und befiehlt, dass das Paket sich in Richtung Mauer bewegt.

Am Rande des Spielfelds wäre es ideal, die Gleichungen so hinzubiegen, dass austretende Wellen sauber rückstandsfrei aufgesaugt werden. Das ist nicht so leicht. Man könnte etwa alle Ränder mit Absorptions- Schichten versehen, mehrere Wellenlängen dick, frei von Reflektionen. (Sowas wie die Wandverkleidung eines schalltoten Raums.) Jedoch die einfachsten Randbedingungen sind entweder, die Wellenfunkton am Rand auf Null zu setzen oder sie stetig periodisch fortzusetzen. Beim Erzwingen der Null werden auslaufende Wellen voll reflektiert. Im Fall der periodischen Randbedingung treten solche Wellen statt dessen an der gegenüber liegenden Kante wieder ein. Das Skript macht es sich einfach und nutzt periodische Randbedingungen.

Wegen des Gauß-Profils von nur wenigen Wellenlängen ist zu erwarten, dass kein scharf definierter Impuls-Eigenwert vorliegt und dass daher die Beugung nach dem Doppelspalt keine völlig dunklen Zonen hergibt, in denen die Amplitude durch destruktive Interferenz ganz auf Null geht. In der Tat. Es sind aber das mittlere Maximum der Amplitude und zwei seitliche Keulen klar zu sehen, also das erwartete Verhalten bei der linearen Überlagerung von Wellen. Das Skript schreibt alle 20 Zeitschritte eine SVG-Datei, in der einige Linien der Niveaus von Wellenamplituden gezeichnet sind. Je dicker und dunkler die Linie, um so größer die Amplitude.

Die ausgegebene Dateien heißen: wave-0.svg ... wave-12.svg

Am Ende der Simulation gibt es noch eine Datei "diffrac.svg", in der die Linien für das zeitgemittelte Amplitudenquadrat allein in der Zone hinter dem Doppelspalt auftauchen. Die drei Keulen der maximalen Intensität treten klar hervor. In der hier gewählten Darstellung liegt der Doppelspalt am unteren Bildrand.

In den ersten Phasen der Simulation prallt das Wellenpaket auf die Mauer und entwickelt ein Knautsch-Muster, worin sich die ankommende und die reflektierte Welle überlagern. Im weiteren Verlauf wandert die große Masse des Pakets nach dem Aufprall zurück und streut auch in andere Richtungen. Der kleine Anteil, der es durch den Doppelspalt schafft, entfaltet die Maxima und Minima von Amplituden, die wie erwartet von einer Interferenz zeugen.

Wie üblich sind die Gleichungen für die Numerik dimensionslos aufbereitet. Anders ausgedrückt, im vorliegenden Beispiel werden das Wirkungsquantum und die Teilchenmasse hinweg skaliert: ℏ=1, (2m)=1. Der Spitzenwert von V auf dem Potenzialwall beträgt 1. Die Schritte d,h für die Raster von Raum und Zeit haben hier auch den überaus praktischen Wert 1.

${\displaystyle i\partial _{t}f=H(f)=-(\partial _{x}^{2}+\partial _{y}^{2})f+v\cdot f}$

${\displaystyle H=H_{x}+H_{y};\quad H_{x}=-\partial _{x}^{2}+(v(x,y)/2);\quad H_{y}=-\partial _{y}^{2}+(v(x,y)/2)}$

Sei der Schritt Δx=Δy=d und Δt=h. Dann gibt es den naiven expliziten numerischen Ansatz, um die Welle bei (t+h) zu bestimmen, wenn sie bei (t) bekannt ist:

${\displaystyle i\cdot (f(t+h)-f(t))/h=(Hf)(t).}$

Oder den besseren, symmetrischen, hier durchgerechneten Ansatz

${\displaystyle i\cdot (f(t+h)-f(t))/h=[(Hf)(t+h)+(Hf)(t)]/2.}$

Diese Methode nach Crank-Nicholson sorgt oft für stabile numerische Iterationen, während mit der naiven Gleichung parasitäre Amplituden sich exponentiell aufschaukeln können.

In den Operatoren ${\displaystyle H_{x},H_{y}}$ werden die zweiten Ableitungen durch zweite Differenzen angenähert nach dem Muster

${\displaystyle \partial _{x}^{2}f(x,y)=(f(x+d,y)+f(x-d,y)-2f(x,y))/d^{2},}$

${\displaystyle \partial _{y}^{2}f(x,y)=(f(x,y+d)+f(x,y-d)-2f(x,y))/d^{2}.}$

Die Wellenfunktion ist eine Tabelle ${\displaystyle F(j,k)=f(x_{j},y_{k})}$ mit zwei Indizes. Sei genauer F die Tabelle für f(t,x,y), G die Tabelle für f(t+h/2,x,y). Der Halbschritt mit dem Operator ${\displaystyle H_{x}}$ bringt dann für jeden festgehaltenen Index k folgende Gleichungen zwischen F,G:

${\displaystyle i\cdot (G(j,k)-F(j,k))/(h/2)={\tfrac {1}{2}}[(G(j+1,k)+G(j-1,k)-2G(j,k))/d^{2}+}$

${\displaystyle +(F(j+1,k)+F(j-1,k)-2F(j,k))/d^{2}+v(j,k)\cdot (G(j,k)+F(j,k))/2]}$

F ist bekannt, so dass G als Funktion des Index j ein tridiagonales lineares Gleichungssystem erfüllt. Bei periodischen Randbedingungen, tridiagonal verziert mit Zahlen unten links und oben rechts in der Matrix. Denn F(N,k) etwa ist mit F(0,k) zu identifizieren (j=0,1,...N-1).

Ein weiterer Halbschritt nach dem gleichen Strickmuster geht von f(t+h/2,x,y) nach f(t+h,x,y) mit dem Operator ${\displaystyle H_{y}}$.

Die wichtigen Teile des Skripts kurz benannt:

Die Klasse class svgdump und die Funktionen myplot2( ), somecontours( ) poduzieren SVG-Grafik-Dateien aus den Tabellen.

Funktion tridiper( ) löst tridiagonale Gleichungen mit periodischer Randbedingung, tridistd( ) ohne solche.

Funktion schroedingerstep( ) macht einen elementaren Schritt bei eindimensionalen und schroedinger2d( ) bei zweidimensionalen Wellen.

Das initiale Gauß-Wellenpaket wird von wavepkt( ) hergestellt und obstaclex( ) baut die Potenzial-Mauer mit zwei Öffnungen.

Hauptfunktion wavetest( ) spielt die ganze Sequenz durch und spuckt die grafischen Darstellungen aus.

Als Zugabe hängt an wavetest2( ) für Geduldige. Lange Laufzeit, etwas "schönere" Bilder. Ein Absorptions-Rahmen für die auslaufenden Wellen umringt das Spielfeld, seine Ausdehnung beträgt zwei Wellenlängen.

# -*- coding: utf-8 -*-
from math import sqrt, pi, sin, cos, exp, sinh, cosh, atan2
#from random import random as rnd # rnd() is in interval [0,1)

def div(p,q) : return int(p/q)
def tanh(x): return sinh(x)/cosh(x)
def coth(x) : return cosh(x)/sinh(x)

def expc(z) : # complex exp
  u=exp(z.real); a=z.imag; c=cos(a); s=sin(a); i=0+1j
  return u*(c+i*s)

def abs2(z) : a,b= z.real,z.imag; return(a*a+b*b)
def phase(z) : # range 0 to 2*pi
  p=atan2(z.imag,z.real); return (p if (p>=0) else (p+2*pi))

def sqrtc(z) :
  a=phase(z)/2; r=sqrt(sqrt(abs2(z))); return r*(cos(a)+1j*sin(a))

def vsum(x,y,f,g) : # Vektorsum, Koeffs f,g
  n=len(x); z=[0]*n
  for i in range(n) : z[i]= f*x[i]+g*y[i]
  return z

def vscale(x,f) : # Vektor-Skalierung mit Faktor
  n=len(x); z=[0]*n
  for i in range(n) : z[i]= f*x[i]
  return z

def aprx(x) : # 3-digit approximation
  y= (-x) if(x<0) else x; z= int(1e3*y+1e-3); z= (-z) if (x<0) else z
  return z/1000.0

def top125(y) :
  z=1e-9; y= max(abs(y),z)
  while z<y : z=10*z
  if (z/5.0)>y : z=z/5.0
  elif (z/2.0)>y : z=z/2.0
  return z

###### primitive plot stuff

class svgdump :
  def plotbegin(sd, w,h,f) : # args: width height zoomfactor
    W=str(w); H=str(h); u=''; fw=str(f*w); fh=str(f*h)
    s=( '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n'+
     '<!DOCTYPE svg>\n'+
     '<svg xmlns="http://www.w3.org/2000/svg"\n'+
     '  xmlns:xlink="http://www.w3.org/1999/xlink"\n'+
     '  width="'+W+'" height="'+H+'" viewBox="0 0 '+fw+' '+fh+'">\n')
    return s

  def plotend(sd) : return '</svg>\n'
  def move(sd,x,y) : return 'M'+str(x)+','+str(y)
  def line(sd,x,y) : return 'L'+str(x)+','+str(y)

  def ball(sd, x,y,r,col) :
    return ('<circle fill="'+col+'" cx="'+str(x)+'" cy="'+str(y)+
      '" r="'+str(r)+'" />\n')
      
  def labl(sd,s,x,y, size=14,fill='#000') :
    f=' font-family="Arial"'
    t='<tspan x="'+str(x)+'" y="'+str(y)+'">'+s+'</tspan>'
    fsz= str(size)
    return ('<text fill="'+fill+'"'+f+' font-size="'+fsz+
      'px">'+t+'</text>\n')

  def path(sd,w,s,lc='#000',fc='none') :
    # path w=with lc=linecolor fc=fillcol
    t='<path fill="'+fc+'" stroke="'+lc+'" stroke-width="'+str(w)+'" '
    t+= ' stroke-linecap="round" ' 
    return t + 'd="\n'+s+'" />\n'

  def curve(sd,scale,xy, fill='#000', width=1, bkgr='none') :
    (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= div(len(xy),2) 
    for i in range(n) : 
      x= aprx(xmin+fx*(xy[2*i]-xs)); y= aprx(ymin+fy*(xy[2*i+1]-ys))
      s+= sd.move(x,y) if(i==0) else sd.line(x,y)
    return sd.path(width,s,lc=fill,fc=bkgr)

  def grids(sd,scale, gx, gy, fill='#000', width=1) :
    (xs,fx,xmin, ys,fy,ymin) = scale; u=''
    (x0,dx,x1, ya,yb) = gx; (y0,dy,y1, xa,xb)= gy # in user units
    ya= aprx(ymin+fy*(ya-ys)); yb= aprx(ymin+fy*(yb-ys))
    xa= aprx(xmin+fx*(xa-xs)); xb= aprx(xmin+fx*(xb-xs))
    while x0<(x1-0.1*dx) : 
      x= aprx(xmin+fx*(x0-xs))
      s= sd.move(x,ya) + sd.line(x,yb); x0+=dx;
      u+= sd.path(width,s,lc=fill,fc='none')
    while y0<(y1-0.1*dy) : 
      y= aprx(ymin+fy*(y0-ys))
      s= sd.move(xa,y) + sd.line(xb,y); y0+=dy
      u+= sd.path(width,s,lc=fill,fc='none')
    return u

  def frame(sd,scale, xmin,ymin, xmax,ymax, dx,dy) :
    xyf=[xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin]
    bf= sd.curve(scale, xyf, fill='#777',bkgr='#ffd',width=3) 
    if (dx>0)and(dy>0)and(xmax>xmin)and(ymax>ymin) :
      gx=(xmin+dx,dx,xmax, ymin,ymax); gy=(ymin+dy,dy,ymax, xmin,xmax)
      bf+= sd.grids(scale, gx, gy, fill='#ccc', width=2) 
    return bf

  def curveset(sd,scale,xyz, fill='#000', width=1, bkgr='none') :
    (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= len(xyz) 
    #print('curveset n='+str(n))
    for i in range(n) : 
      px,py,pz= tuple(xyz[i])
      x= aprx(xmin+fx*(px-xs)); y= aprx(ymin+fy*(py-ys))
      s+= sd.move(x,y) if(pz==0) else sd.line(x,y)
    return sd.path(width,s,lc=fill,fc=bkgr)

  def ballset(sd,scale,plt, fill='#000', r=5) :
    (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= len(plt)
    for i in range(n) : 
      x= aprx(xmin+fx*(plt[i][0]-xs)); y= aprx(ymin+fy*(plt[i][1]-ys))
      s+= sd.ball(x,y,r,fill)
    return s

  def ballsrf(sd,scale,plt) : # vary radius and fillcolor
    (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= len(plt)
    for i in range(n) : 
      x= aprx(xmin+fx*(plt[i][0]-xs)); y= aprx(ymin+fy*(plt[i][1]-ys))
      s+= sd.ball(x,y,plt[i][2],plt[i][3])
    return s

  def label(sd,scale,txt,tx,ty) : # tx,ty user coordinates not pixels
    (xs,fx,xmin, ys,fy,ymin) = scale 
    x= aprx(xmin+fx*(tx-xs)); y= aprx(ymin+fy*(ty-ys))
    return sd.labl(txt,x,y)
  
  def infobox(sd,s,x,y) :
    t=s.split('\n'); k=len(t); u=''
    for i in range(k) : u+= sd.labl(t[i],x,y+i*26,size=24)
    return u

# end class svgdump

def myplot2(fn, xa,ya,xb,yb,dx,dy, pball, pcurve, s, style=0) :
  # dots(x,y,r,color) in pball array. window size 3:2 
  winx=1000;winy=667; dwx=20;dwy=20 # winx=1000;winy=1000;dw=25
  xmin=xa;ymin=ya;xmax=xb;ymax=yb 
  xyf=[xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin]
  xs,fx,xpix= xmin, (winx-2*dwx)/(xmax-xmin+0.0),dwx
  # scale user units, xpix,ypix is pixel-pos at user-min.
  ys,fy,ypix= ymin,-(winy-2*dwy)/(ymax-ymin+0.0),winy-dwy
  scale=(xs,fx,xpix, ys,fy,ypix)
  colo=['#777','#000','#00a','#077','#0a0','#770','#a00']
  wid= [3]*len(colo); ch='fedccba98765'
  if (style==1)or(style==2) : # adapt colors,widths to pcurve 
    n=len(pcurve)+1; colo=['']*n; wid=[1]*n
    for k in range(n): wid[k]=2*k+1; colo[k]='#3'+ch[k]+ch[k]
  if style==2 : colo[1]='#000'; wid[1]=2 # barrier picture
  sd=svgdump(); buf=sd.plotbegin(winx,winy,1)
  buf+= sd.frame(scale, xa,ya, xb,yb, dx,dy)
  buf+= sd.ballsrf(scale,pball)
  for k in range(len(pcurve)) : 
    j=(k+1) % (len(colo))
    buf+= sd.curveset(scale,pcurve[k], fill=colo[j], width=wid[j])
  if len(s)>0 : buf+= sd.infobox(s,2*winx/3,2*winy/3)
  g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()

def contourset(z,a) : # level contours z(x,y)=a
  #
  def closecontours(p) :
    n=len(p); q=[]; ix=-1; n1=n-1
    for i in range(n) :
      stop=(i==n1)or(p[i+1][2]==0); start=(p[i][2]==0); q+=[p[i]]
      if stop and (ix>=0)and(ix!=i) :
        d=abs(p[i][0]-p[ix][0])+abs(p[i][1]-p[ix][1])
        if d<2 : q+= [ [p[ix][0],p[ix][1],1] ]
      if start: ix=i
    return q  
  # 
  nx=len(z); ny=len(z[0]); nx1=nx-1; ny1=ny-1; pts=[]; a=0.0+a #float
  for j in range(nx1) :
    for k in range(ny1) : 
      p=z[j][k]-a; q=z[j][k+1]-a; r=z[j+1][k]-a
      if (p*q)<0 : t=p/(p-q); pts+=[[j,k+t,0]]
      if (p*r)<0 : t=p/(p-r); pts+=[[j+t,k,0]]
  n=len(pts); m=0; j=0; z=0; plt=[[]]*n; pts[j][2]=1
  print('contourset a='+str(a)+' n='+str(n))
  while m<n : # nb of plot points
    d=nx+ny; x=pts[j][0]; y=pts[j][1]; ix=-1; plt[m]=[x,y,z]; m+=1
    for k in range(n) :
      dk=abs(pts[k][0]-x)+abs(pts[k][1]-y) 
      if (dk<d) and (pts[k][2]==0) :  d=dk; ix=k
    j=ix; z=1 if (d<2) else 0; pts[j][2]=1
  return closecontours(plt) # list of x,y,z. z=0 pen up, z=1 pen down.
    
def somecontours(z,m) : # make m level contours for amplitude array z
  nx=len(z); mx=0
  for j in range(nx) : u=max(z[j]); mx=max(mx,u); cs=[[]]*m
  for k in range(1,m+1) : cs[k-1]=contourset(z, k*mx/(m+1.0))
  return cs # curve set to be used in myplot2()

### end plotting
 
def tridicheck(x,a,b,c,y, dbg=0) : # trdiagonal equations
  n=len(x); n1=n-1; dmax=0; imax=0; yq=0; z=0; dv=[0]*n
  for i in range(n) : 
    yp= 0 if(i>=n1) else y[i+1]
    u=x[i]-a[i]*yq-b[i]*y[i]-c[i]*yp; d=abs2(u);z+=d;yq=y[i];dv[i]=-u
    if d>dmax : dmax=d; imax=i
  if dbg>1 : print('tridiag dbg='+str(dbg)+
    ' err='+str(z)+' imax='+str(imax)+'/'+str(n))
  return z,dv # error and difference vector

def tridistd(x,a,b,c) : # xi=ai*y(i-1)+bi*yi+ci*y(i+1)
  n=len(x); bb=[0]*n; xx=[0]*n; y=[0]*n; bb[0]=b[0]; xx[0]=x[0]; z=0
  for i in range(1,n) :
    f=a[i]/bb[i-1]; bb[i]=b[i]-f*c[i-1]; xx[i]=x[i]-f*xx[i-1] 
  n1=n-1; y[n1]= xx[n1]/bb[n1]; i=n1; h=n1-1; yp=y[n1] 
  while h>=0 : y[h]=(xx[h]-c[h]*yp)/bb[h]; yp=y[h]; h-=1
  return y

def tridiper(x,a,b,c) : # tridiagonal-periodic-boundary equation
  #
  def solve2(r,s) :
    a,b,u= tuple(r); c,d,v=tuple(s); # ax+by=u; cx+dy=v
    q=a*d-b*c; x=(d*u-b*v)/q; y=(-c*u+a*v)/q; return x,y 
  #
  n=len(x); y=[0]*n; b1=[0]*n;x1=[0]*n; b1[0]=b[0];x1[0]=x[0]
  a1=[0]*n; a1[0]=a[0] # last column gets populated
  for i in range(1,n) : # get (b1[n1]+a1)*y[n1]+c[n1]*y[0]=x1[n1]
    f=a[i]/b1[i-1]; b1[i]=b[i]-f*c[i-1];
    x1[i]=x[i]-f*x1[i-1]; a1[i]=-f*a1[i-1] 
  b2=[0]*n; x2=[0]*n; n1=n-1;i=n1-1; b2[n1]=b[n1];x2[n1]=x[n1];c2=c[n1]
  while i>=0 : # get (b2[0]+c2)*y[0]+a[0]*y[n1]=x2[0]
    f=c[i]/b2[i+1]; b2[i]=b[i]-f*a[i+1];
    x2[i]=x[i]-f*x2[i+1]; c2=-f*c2; i-=1
  y[0],y[n1]=solve2([c[n1],b1[n1]+a1[n1],x1[n1]],[b2[0]+c2,a[0],x2[0]])
  h=n1-1; yp=y[n1]; yq=yp; z=0;
  while h>=1 : y[h]=(x1[h]-c[h]*yp-a1[h]*yq)/b1[h]; yp=y[h]; h-=1
  for i in range(n) : 
    yp= y[0] if(i>=n1) else y[i+1]
    u=x[i]-a[i]*yq-b[i]*y[i]-c[i]*yp; d=abs2(u); z+=d; yq=y[i]
  return y,z

''' schroedinger equ Crank-Nicholson style
2D Schroedinger plane waves
  i dot(f)=Hf= -Delta f + V*f = [V[i]*f[i]+(2f[i]-f(i+1]+f[i+1])/hh]
  i (g-f)/t = (H g + H f)/2 ;  t=timestep, h=space step.
  g-f=(-it/2)(Hg+Hf); (1+itH/2)g=(1-itH/2)f to be solved for g.
  alternate directions, operator split.
'''

def schroedingerstep(f,v,dt,dx,fv,bc=1) : # transform f to new array g
  #bc=1 # periodic boundary
  #bc=0 # zero boundary cond
  n=len(f); h2=dx*dx; i=0+1j # solve (1+itH/2)g=(1-itH/2)f
  a=[0]*n; b=[0]*n; c=[0]*n; d=[0]*n; z=i*dt/2.0; za=0; zb=0
  for j in range(n) :
    ff=f[j]; fm=f[j-1] if (j>0) else za; fp=f[j+1] if (j<(n-1)) else zb
    if bc==1 :
      jm=(j-1) if (j>0) else (n-1); jp=(j+1) if (j<(n-1)) else 0
      d[j]= ff-z*(v[j]*fv*ff+(2*ff-f[jp]-f[jm])/h2)
    else : d[j]= ff-z*(v[j]*fv*ff+(2*ff-fp-fm)/h2)
    a[j]= -z/h2; c[j]= -z/h2
    b[j]= 1+z*(v[j]*fv+2.0/h2)
  if bc==1 : g,err= tridiper(d,a,b,c)
  else : g= tridistd(d,a,b,c); err,dv=tridicheck(d,a,b,c,g)
  return g,err

def schroedinger2d(f,v,dt,dx,bc=1) : # operator split, 2 halfsteps
  nx=len(f); ny=len(f[0]); g=[[]]*nx; es=0; fv=0.5 # potential split
  for ix in range(nx) :
    g[ix],err = schroedingerstep(f[ix],v[ix],dt/2.0,dx,fv,bc); es+=err
    if err>0.1 : print('s2d err='+str(err)+' ix='+str(ix))
  ff=[0]*nx; vv=[0]*nx; gg=[0]*nx
  for iy in range(ny) :
    for ix in range(nx) : ff[ix]=g[ix][iy]; vv[ix]=v[ix][iy]
    gg,err = schroedingerstep(ff,vv,dt/2.0,dx,fv,bc); es+=err
    if err>0.1 : print('s2d err='+str(err)+' iy='+str(iy))
    for ix in range(nx) : f[ix][iy]=gg[ix]
  return es # es=errorsum; array f is modified in place

def wavepkt(nx,ny, cx,cy, dx,dy, kx,ky) : # define wavepacket f
  f=[[]]*nx; i=0+1j
  for ix in range(nx) :
    f[ix]=[0+0j]*ny
    for iy in range(ny) :
      x=(ix-cx)/(dx+0.0); y=(iy-cy)/(dy+0.0)
      f[ix][iy]= exp(-x*x-y*y) * expc(i*(ix*kx+iy*ky))
  return f

def obstaclex(nx,ny, vv,px,dy) : # double slit at px, y-separation dy
  v=[[]]*nx; ya=div(ny-dy,2); yb=div(ny+dy,2)
  wall=[vv/3.0,2*vv/3.0,vv,2*vv/3.0,vv/3.0] # 5 pixels
  for ix in range(nx) : v[ix]=[0]*ny
  for j in range(5) :
    qx=px-2+j
    for iy in range(ny) :
      ok= (iy<(ya-2))or(iy>(yb+2)) or ((iy>(ya+2))and(iy<(yb-2)))
      if ok or (dy<0) : v[qx][iy]= wall[j]
  return v

def wavetest( nx=120, ny=90, nt=240) : # double slit diffraction
  #
  def waveplot(fn,f,v, mode=0) : # wave amplitude contours
    fnm=0; mx=0; col='#7f0'; clv='#07f'; nx=len(f); ny=len(f[0])
    balls=[]; xm=0;ym=0; gx=0;gy=0; vm=0; amp=[[]]*nx
    for ix in range(nx) : amp[ix]=[0]*ny
    for ix in range(nx) :
      for iy in range(ny) :
        z=abs2(f[ix][iy]); fnm+=z; gx+=z*ix;gy+=z*iy
        vm=max(vm, v[ix][iy]); amp[ix][iy]=sqrt(z)
        if z>mx : mx=z; xm=ix; ym=iy
    print('Maxval='+str(mx)+'@'+str([xm,ym])+' Norm2='+str(fnm))
    print('Center @'+str([gx/fnm,gy/fnm])) # center of gravity
    cv= [] if (vm==0) else somecontours(v,1)
    cs=somecontours(amp,5); curves= cv+cs if (mode==1) else []
    for ix in range(nx) :
      for iy in range(ny) :
        if vm>0 : r=int(50.0*v[ix][iy]/vm)/10.0
        if (mode==0) and (r>=1) : balls+=[[ix,iy, r,clv]] 
        r=sqrt(abs2(f[ix][iy])/mx)*6; r=int(10*r)/10.0
        if (mode==0) and (r>=1) : balls+=[[ix,iy, r,col]] 
    if mode==1 : balls=[]
    else : print('len(balls)='+str(len(balls)))  
    myplot2(fn, 0,0, nx,ny, 10,10, balls,curves,'',style=2)
  #
  def ampsquareadd(a,f,m) : # amplitudes after the slits
    n=len(f); ny=len(f[0])
    for i in range(m,n) :
      if len(a[i])==0 : a[i]=[0]*ny # first pass
      for k in range(ny): a[i][k]+= abs2(f[i][k])
  #
  def diffracplot(fn,a) : # summed-up diffraction pattern
    na=len(a); nb=len(a[na-1]); m=na-1
    while len(a[m])>0 : m-=1
    m+=1; ny=na-m; nx=nb; p=[[]]*nx; mx=0
    for i in range(nx) :
      p[i]=[0]*ny
      for k in range(ny) : z=a[k+m][i]; p[i][k]=sqrt(z); mx=max(mx,z)
    print('ampsquare max='+str(mx))
    cs=somecontours(p, 8)
    myplot2(fn, 0,0, nx,ny, 10,10, [],cs,'',style=1)
  #
  cx=div(nx,3);cy=div(ny,2); dx=div(nx,4);dy=div(ny,3); lambd=20 
  kx=2*pi/lambd; ky=0; f=wavepkt(nx,ny, cx,cy, dx,dy, kx, ky)
  mode=1; dk=20; vmax=1 #  vmax=1 potential wall height
  dsy=div(3*lambd,2) # double slit distance
  v=obstaclex(nx,ny,vmax, 2*cx, dsy); dx=1; dt=1; asq=[[]]*nx
  for k in range(nt+1) : 
    # packet move ~= 0.3 units in positive x direction.
    es=schroedinger2d(f,v,dt,dx) 
    print('t='+str(k)+' errsum='+str(es))
    if (k%dk)==0 :  waveplot('wave-'+str(k/dk), f,v,mode)
    if k> div(nt,3) : ampsquareadd(asq,f,2*cx+10) 
  diffracplot('diffrac', asq)

wavetest( nx=120, ny=90, nt=240)

def wavetest2( nx=300,ny=200, nt=600, lambd=20, imi=10) :
  #wavelength 20 dots, image interval 20 ticks
  #
  def waveplot(fn,f,v,d) : # wave amplitude contours, inside frame d
    fnm=0; mx=0; nx=len(f); ny=len(f[0]); nxd=nx-2*d; nyd=ny-2*d
    vv=[[]]*nxd; xm=0;ym=0; gx=0;gy=0; vm=0; amp=[[]]*nxd
    for ix in range(nxd) : amp[ix]=[0]*nyd; vv[ix]=[0]*nyd
    for ix in range(nxd) :
      for iy in range(nyd) :
        z=abs2(f[ix+d][iy+d]); fnm+=z; gx+=z*(ix+d);gy+=z*(iy+d)
        u=v[ix+d][iy+d];vm=max(vm,u); amp[ix][iy]=sqrt(z);vv[ix][iy]=u
        if z>mx : mx=z; xm=ix; ym=iy
    print('Maxval='+str(mx)+'@'+str([xm,ym])+' Norm2='+str(fnm))
    print('Center @'+str([gx/fnm,gy/fnm])) # center of gravity
    cv= [] if (vm==0) else somecontours(vv,1)
    cs=somecontours(amp,5); curves= cv+cs
    myplot2(fn, 0,0, nx-2*d,ny-2*d, 10,10, [], curves,'',style=2)
  #
  def ampsquareadd(a,f,m) : # amplitudes after the slits
    n=len(f); ny=len(f[0])
    for i in range(m,n) :
      if len(a[i])==0 : a[i]=[0]*ny # first pass
      for k in range(ny): a[i][k]+= abs2(f[i][k])
  #
  def diffracplot(fn, a, d) : # diffraction pattern in frame(d)
    na=len(a); nb=len(a[na-1]); m=na-1
    while len(a[m])>0 : m-=1
    m+=1; ny=na-m-d; nx=nb; p=[[]]*(nx-2*d); mx=0
    for i in range(nx-2*d) :
      p[i]=[0]*ny
      for k in range(ny) : z=a[k+m][i+d]; p[i][k]=sqrt(z); mx=max(mx,z)
    print('ampsquare max='+str(mx))
    cs=somecontours(p, 8)
    myplot2(fn, 0,0, nx-2*d,ny, 10,10, [],cs,'',style=1)
  #
  def absorber(f,n) : # use with zero boundary cond. v[0]=f,v[n-1]=1
    v=[0]*n; n1=n-1.0; g=1-f
    for i in range(n) : z=(cos(i*pi/n1)+1)/2; v[i]=1-g*z
    return v
  #
  def absorb(f,a) : # after each iteration, absorber frame damps wave
    nx=len(f); ny=len(f[0]); m=len(a)
    for k in range(nx) : 
      for j in range(m) : f[k][ny-1-j] *= a[j]; f[k][j] *= a[j]
    for k in range(ny) : 
      for j in range(m) : f[nx-1-j][k] *= a[j]; f[j][k] *= a[j]
  #
  cx=div(nx,4);cy=div(ny,2); dx=div(nx,6);dy=div(ny,4) #lambd=20 
  kx=2*pi/lambd; ky=0; f=wavepkt(nx,ny, cx,cy, dx,dy, kx,ky)
  dk=imi; vmax=1 #  dk=plot interval, vmax=1 potential wall height
  dsy= 2*lambd  #  double slit distance
  v=obstaclex(nx,ny,vmax, 2*cx, dsy); dx=1; dt=1; asq=[[]]*nx
  nm=2*lambd; ab=absorber(0.9,nm) # 0.9=damping, nm=absorber margins
  for k in range(nt+1) : 
    # packet move ~= 0.3 units per step in positive x direction.
    es=schroedinger2d(f,v,dt,dx, bc=0) 
    print('t='+str(k)+' errsum='+str(es))
    absorb(f,ab)
    if (k % dk)==0 :  waveplot('wave2-'+str(k/dk), f,v,nm)
    if k> div(nt,3) : ampsquareadd(asq,f,2*cx+10) 
  diffracplot('diffrac2',asq, nm)

wavetest2()

Um das Ergebnis animiert abzuspielen, kann ein Shellskript folgender Art dienen (ist hier abhängig von: rsvg-convert,ffmpeg,ffplay). Die 60 Bilder werden ins Format .png umgemodelt und dann im Takt von einem pro Sekunde zu einem Mini-Video verschmolzen.

# animation, matroska video
if [ ! -f waves.mkv ]; then
  for X in wave2-*.svg; do Y=${X%.svg}; rsvg-convert $X >$Y.png ; done
  ffmpeg -f image2 -r 1 -i wave2-%d.png -s vga -r 1 waves.mkv
fi
ffplay waves.mkv

Hoch zum Seitenanfang