harmoszop.html

Ein Wellenpaket in die SGL einsetzen...

`>`	restart:

`>`	*sgl:=diff(u(x),x$2)=-2(E-V(x))u(x);*

sgl := diff(u(x),`$`(x,2)) = -2*(E-V(x))*u(x)

Wellenpaket (Normierung sqrt(Pi)):

`>`	gp:=exp(-x^2/2);

gp := exp(-1/2*x^2)

`>`	subs(u(x)=gp,sgl);

diff(exp(-1/2*x^2),`$`(x,2)) = -2*(E-V(x))*exp(-1/2*x^2)

`>`	simplify(%);

exp(-1/2*x^2)*(-1+x^2) = -2*(E-V(x))*exp(-1/2*x^2)

`>`	expand(%/2/gp);

-1/2+1/2*x^2 = -E+V(x)

Für ein quadratisches Potential ist also das Gaußpaket zur Energie E = 1/2 eine Lösung der stationären SGL. Man kann das Wellenpaket "anhalten"!

`>`	V:=x->x^2/2;

Operatoren

Für den harmonischen Oszillator (quadratisches Potential) lässt sich die SGL auch so schreiben:

`>`	*(x^2u(x)-diff(u(x),x$2))/2-Eu(x)=0;*

1/2*x^2*u(x)-1/2*diff(u(x),`$`(x,2))-E*u(x) = 0

Das erinnert die dritte binomische Formel - hier allerdings für Operatoren

`>`	dop:=Diff(``,x): dop2:=Diff(``,x$2):

>

`>`	*x^2-dop2=(x+dop)(x-dop);**

x^2-Diff(``,`$`(x,2)) = (x+Diff(``,x))*(x-Diff(``,x))

>

`>`	*Op:=u->((xu+diff(u,x))/sqrt(2));**

Op := proc (u) options operator, arrow; (x*u+diff(u,x))/sqrt(2) end proc

`>`	*Ops:=u->(xu-diff(u,x))/sqrt(2);**

Ops := proc (u) options operator, arrow; (x*u-diff(u,x))/sqrt(2) end proc

>

Anwendung des Operators

`>`	simplify((Ops@Op)(u(x)));

1/2*x^2*u(x)-1/2*u(x)-1/2*diff(u(x),`$`(x,2))

Das ist die linke Seite der Schrödingergleichung zur Energie 1/2.

`>`	simplify((Op@Ops)(u(x)));

1/2*x^2*u(x)+1/2*u(x)-1/2*diff(u(x),`$`(x,2))

Das ist etwas anderes, es kommt also auf die Reihenfolge der Anwendung der Operatoren an und es gilt:

> %-%%;

Der Kommutator der beiden Operatoren ist also der Einheitsoperator

>

Was machen die Operatoren aus dem Paket?

`>`	(Op@Ops)(gp);

exp(-1/2*x^2)

Der Grundzustand bleibt erhalten

`>`	(Ops@Op)(gp);

oder wird vernichtet. Wie funktioniert das?

`>`	Ops(gp);

x*exp(-1/2*x^2)*2^(1/2)

Eine neue Funktion? Und wieder zurück:

`>`	Op(%);

exp(-1/2*x^2)

>

`>`	subs(u(x)=Ops(gp),sgl);

diff(x*exp(-1/2*x^2)*2^(1/2),`$`(x,2)) = -2*(E-1/2*x^2)*x*exp(-1/2*x^2)*2^(1/2)

`>`	simplify(%);

x*exp(-1/2*x^2)*2^(1/2)*(-3+x^2) = (-2*E+x^2)*x*exp(-1/2*x^2)*2^(1/2)

`>`	expand(%/gp/sqrt(2));

>

Das stimmt für E = 3/2! Ops erzeugt also aus dem Grundzustand eine neue Lösung. Wir vergleichen mit den bekannten Lösungen (ohne sqrt(Pi)):

>

`>`	*v:=(n,x)->1/sqrt(2^nn!)HermiteH(n,x)exp(-x^2/2); #Pi^(-1/4);*

v := proc (n, x) options operator, arrow; 1/sqrt(2^n*n!)*HermiteH(n,x)*exp(-1/2*x^2) end proc

`>`	simplify(v(1,x));

x*exp(-1/2*x^2)*2^(1/2)

`>`	Ops(gp);

x*exp(-1/2*x^2)*2^(1/2)

>

`>`	simplify((Ops@@3)(gp));

2^(1/2)*x*exp(-1/2*x^2)*(2*x^2-3)

`>`	simplify(v(3,x));

1/3*3^(1/2)*x*(2*x^2-3)*exp(-1/2*x^2)

`>`	#int(gp^2,x=-infinity..infinity);

`>`	seq(int((Ops@@n)(gp)^2,x=-infinity..infinity),n=1..5);

Zur Normierung fehlt also nur der Faktor

`>`	1/sqrt(n!)/Pi^(1/4);

1/(n!^(1/2)*Pi^(1/4))

Wir verzichten wieder auf sqrt(Pi)

>

`>`	erz:=(u,n)->(Ops@@n)(u)/sqrt(n!);

erz := proc (u, n) options operator, arrow; `@@`(Ops,n)(u)/sqrt(n!) end proc

`>`	test:=simplify(erz(gp,6));

test := 1/60*exp(-1/2*x^2)*(8*x^6-60*x^4+90*x^2-15)*5^(1/2)

>

`>`	simplify(v(6,x));

1/60*exp(-1/2*x^2)*(8*x^6-60*x^4+90*x^2-15)*5^(1/2)

>

`>`	vern:=(u,n)->(Op@@n)(u)/sqrt(n!);

vern := proc (u, n) options operator, arrow; `@@`(Op,n)(u)/sqrt(n!) end proc

`>`	simplify(vern(test,6));

exp(-1/2*x^2)

>

Leiteroperatoren mit Normierung

`>`	auf:=(u,n)->Ops(u)/sqrt(n+1);

`>`	auf(test,6);

1/14*(1/30*x*exp(-1/2*x^2)*(8*x^6-60*x^4+90*x^2-15)*5^(1/2)-1/60*exp(-1/2*x^2)*(48*x^5-240*x^3+180*x)*5^(1/2))*2^(1/2)*7^(1/2)

`>`	simplify(%);

1/420*exp(-1/2*x^2)*5^(1/2)*x*(8*x^6-84*x^4+210*x^2-105)*2^(1/2)*7^(1/2)

>

`>`	ab:=(u,n)->Op(u)/sqrt(n);

`>`	simplify(ab(test,6));

1/30*exp(-1/2*x^2)*x*(4*x^4-20*x^2+15)*5^(1/2)*3^(1/2)

`>`	#simplify(v(5,x));

>

Reihenansatz

`>`	restart:

`>`	with(powseries):

`>`	*sgl:=diff(psi(x),x$2)+2(E-x^2/2)psi(x)=0;*

sgl := diff(psi(x),`$`(x,2))+2*(E-1/2*x^2)*psi(x) = 0

Mit powsolve(DGL)

`>`	sol:=powsolve(sgl);

wird die DGL mit einem Reihenansatz (formale Reihe) gelöst:

`>`	tpsform(sol,x,7);

series(C0+C1*x+(-E*C0)*x^2+(-1/3*E*C1)*x^3+(1/6*E^2*C0+1/12*C0)*x^4+(1/30*E^2*C1+1/20*C1)*x^5+(-1/15*E*(1/6*E^2*C0+1/12*C0)-1/30*E*C0)*x^6+O(x^7),x,7)

C0 und C1 sind die Werte für Anfangsbedingungen psi(0) und D(psi)(0).

>

Koeffizienten a(_k) der Reihe

`>`	sol(_k);

Die Reihe sollte abbrechen. Für welches E?

`>`	solve(sol(_k),E);

Die Rekursion stört. Wir versuchen es mit einer Reihe, die schon das GP enthält:

>

`>`	*psi:=x->exp(-x^2/2)f(x);**

psi := proc (x) options operator, arrow; exp(-1/2*x^2)*f(x) end proc

`>`	powsolve(sgl) ;

Error, (in spcldif) final value in for loop must be numeric or character

> sgl;

-exp(-1/2*x^2)*f(x)+x^2*exp(-1/2*x^2)*f(x)-2*x*exp(-1/2*x^2)*diff(f(x),x)+exp(-1/2*x^2)*diff(f(x),`$`(x,2))+2*(E-1/2*x^2)*exp(-1/2*x^2)*f(x) = 0

`>`	# muss erst vereinfacht werden

>

DGL für f(x)

`>`	*sgls:=simplify(simplify(factor(sgl))exp(x^2/2));**

sgls := -f(x)-2*x*diff(f(x),x)+diff(f(x),`$`(x,2))+2*E*f(x) = 0

`>`	sols := powsolve(sgls) :

`>`	sols(_k);

Verschiebung des Index

`>`	subs(_k=n+2,%);

Eigenwerte

`>`	solve(%,E);

`>`	f7:=tpsform(sols,x,8);

`>`	subs(E=7+1/2,f7);

`>`	*simplify(HermiteH(7,x)3/7!);**

Scheint zu passen...

`>`	*#v:=(n,x)->1/sqrt(2^nn!)HermiteH(n,x)exp(-x^2/2); #Pi^(-1/4);*

`>`	#simplify(v(7,x));

Und hier ist die Prozedur:

`>`	restart:

> psi:=proc(n,x)
local sgls,sols,f,u,N;
sgls := -f(x)-2*diff(f(x),x)*x+diff(f(x),`$`(x,2))+2*f(x)*E = 0;
sols := powseries[powsolve](sgls) :
f:=convert(powseries[tpsform](sols,x,n+1),polynom);
f:=subs(C0=Re(I^n),C1=Im(I^n),E=n+1/2,f);
u:=exp(-x^2/2)*f;
N:=int(u^2,x=-infinity..infinity);
u/sqrt(N);
end;

>

`>`	psi(7,x);

1/4*exp(-1/2*x^2)*(-x+2*x^3-4/5*x^5+8/105*x^7)*70^(1/2)/Pi^(1/4)

`>`	*#v:=(n,x)->1/sqrt(2^nn!)HermiteH(n,x)exp(-x^2/2)Pi^(-1/4):*

`>`	#simplify(v(7,x));

>

`>`	plot(psi(7,x),x=-6..6);

[Maple Plot]

>

`>`	plot(psi(8,x)^2,x=-6..6);

[Maple Plot]

Zurück zum Inhalt

komma@oe.uni-tuebingen.de