Prozeduren zur numerischen Lösung der stationären SGL
> | restart:with(plots): |
> | sglp:=proc(V,E,x0,psi0,dx0,dpsi0) local sgl,sol; |
> | sgl:=diff(psi(x),x$2)=-2*(E-V)*psi(x); |
> | sol:=proc(V,E,x0,psi0,dx0,dpsi0) dsolve({sgl,psi(x0)=psi0,D(psi)(dx0)=dpsi0},numeric,output=listprocedure); end proc; rhs(sol(V,E,x0,psi0,dx0,dpsi0)[2]); |
> | end proc; |
> |
Aufruf mit sglp(Potentialterm, Energie, x0, psi(x0), x1, D(psi)(x1))(x)
> | sglp(x^2/2,1,0,1,0,0.1)(1); |
Lineares Potential
> | plot([5*sglp(x,9.473622,0,0,0,1)(x)+9.4735,x,9.4735],x=0..15, color=[red,blue,black],scaling=constrained,thickness=2); |
> |
Quadratisches Potential
> | plot([5*sglp(x^2/2,7/2,0,0,0,1)(x)/4+7/2,x^2/2,7/2],x=-5..5,0..5,color=[red,blue,black],thickness=2); |
> |
> | plot([seq(sglp(x^2/2,E,0,0,0,1)(x),E=seq(3/2+2*n,n=0..5))],x=-6..6,-1..1); |
> |
Coulombpotential ## läuft wegen Singularität ohne den Zusatz 1e-20 nur nach restart (und nicht, wenn vorher andere Potentiale gerechnet wurden)?
> | V:=-1/sqrt(x^2+1e-20): E:=-1/2/9: |
> | plot([sglp(V,E,0,0,0,1)(x)/20+E,V,E],x=-50..50,-0.2..0, color=[red,blue,black],thickness=2,tickmarks=[2,2]); |
> |
Kontinuum (Coulombwellen)
> | V:=-1/sqrt(x^2+1e-20): E:=1/100: # (zumindest) im Kontinuum Randbed. kritisch |
> | plot([sglp(V,E,0,0,0,1)(x)/4+E,V,E],x=-80..80,-1/2..1/4, color=[red,blue,black]); |
> |
Topf
> | V:=Heaviside(-x-1)+Heaviside(x-1); |
> | plot([sglp(200*V-200,-198.8811,0,1,0,0)(x),V],x=-1.4..1.4); |
> |
Potentialstufe
> | V:=x->piecewise(x<10,0,10); |
> | E:=5:x0:=2.300692: |
> | plot([5*sglp(V(x),E,x0,0,x0,1)(x)+E,V(x),E],x=0..15,0..12,numpoints=500,color=[red,blue,black]); |
> |
Wand
> | V:=x->piecewise(x<10,0,x<10.5,10,4); |
> | E:=5:x0:=2.300692: |
> | plot([5*sglp(V(x),E,x0,0,x0,1)(x)+E,V(x),E],x=0..18,0..12,numpoints=500,color=[red,blue,black],thickness=2); |
> |
Kontinuum
> | V:=x->piecewise(x<10,0,x<14,10,4); |
> | E:=12:x0:=0: |
> | plot([5*sglp(V(x),E,x0,0,x0,1)(x)+E,V(x),E],x=0..18,0..15,numpoints=500,color=[red,blue,black]); |
Generelles Problem bei numerischen Verfahren: Randwerte können nicht mehrfach gesetzt werden (z.B. beim Übergang in mehrere Medien).
> |
Rampe
> | V:=x->piecewise(x<0,0,x<10,x,x>10,10); |
> | E:=9.066244: |
> | plot([5*sglp(V(x),E,-10,0,-10,1)(x)+E,V(x),E],x=-5..15,-2..15,numpoints=500,color=[red,blue,black]); |
> |
Keil
> | V:=x->piecewise(x<0,0,x<10,x,x<10.02,10,0); |
> | E:=9.066244: |
> | plot([5*sglp(V(x),E,-10,0,-10,1)(x)+E,V(x),E],x=-5..20,-2..15,numpoints=500,color=[red,blue,black]); |
> |
Periodisches Potential -> Bänder
> | V:=x->signum(sin(a*x)); |
Mögliche Werte für E und Potentialtiefe b (halbe Periode = ganzzahliges Vielfaches von lambda/2):
> | n:='n':m:='m':E:=(n^2+m^2)/4;b:=(m^2-n^2)/4; |
m, n auch halbzahlig
> | a:=1/2:n:=2:m:=9/2:E;b; |
> | plot([sglp(b*V(x),E,0,0,0,E/2)(x)+E,b*V(x),E],x=-8..15,-abs(b)-1/2..E+E/3,numpoints=500, color=[red,blue,black]); |
> |
Bei allen EWn ist "erstaunlich" wie exakt sie angegeben werden müssen, also mit numerischen Methoden unerreichbar sind (sie divergieren immer):
3D-Darstellung über E und x für Oszillator
> |
komma@oe.uni-tuebingen.de