Zurück zum Inhalt

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;

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},nume...
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},nume...
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},nume...
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},nume...
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},nume...
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},nume...
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},nume...
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},nume...

 

>   

 

Aufruf mit sglp(Potentialterm, Energie, x0, psi(x0), x1, D(psi)(x1))(x)

 

>    sglp(x^2/2,1,0,1,0,0.1)(1);

.279555282651855196

 

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);

[Maple Plot]

 

>

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);

[Maple Plot]

 

>   

 

>    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);

[Maple Plot]

 

>   

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]);

[Maple Plot]

 

>   

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]);

[Maple Plot]

 

 

>   

Topf

 

>    V:=Heaviside(-x-1)+Heaviside(x-1);

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);

[Maple Plot]

 

>   

Potentialstufe

 

>    V:=x->piecewise(x<10,0,10);

V := proc (x) options operator, arrow; piecewise(x < 10,0,10) end proc

 

>    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]);

[Maple Plot]

 

>   

Wand

 

>    V:=x->piecewise(x<10,0,x<10.5,10,4);

V := proc (x) options operator, arrow; piecewise(x < 10,0,x < 10.5,10,4) end proc

 

>    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);

[Maple Plot]

 

>   

Kontinuum 

 

>    V:=x->piecewise(x<10,0,x<14,10,4);

V := proc (x) options operator, arrow; piecewise(x < 10,0,x < 14,10,4) end proc

 

>    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]);

[Maple Plot]

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);

V := proc (x) options operator, arrow; piecewise(x < 0,0,x < 10,x,10 < x,10) end proc

 

>    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]);

[Maple Plot]

 

>   

Keil

 

>    V:=x->piecewise(x<0,0,x<10,x,x<10.02,10,0);

V := proc (x) options operator, arrow; piecewise(x < 0,0,x < 10,x,x < 10.02,10,0) end proc

 

>    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]);

[Maple Plot]

 

>   

 Periodisches Potential  -> Bänder

 

>    V:=x->signum(sin(a*x));

V := proc (x) options operator, arrow; signum(sin(a*x)) end proc

 

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;

E := 1/4*n^2+1/4*m^2

b := 1/4*m^2-1/4*n^2

m, n auch halbzahlig

 

>    a:=1/2:n:=2:m:=9/2:E;b;

97/16

65/16

 

>    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]);

[Maple Plot]

 

>   

 

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

 

>   

[Maple Plot]

 

Zurück zum Inhalt

komma@oe.uni-tuebingen.de