hydrogen_10.mws

Moderne Physik mit Maple

PDF-Buch Moderne Physik mit Maple

Update auf Maple 10

Kapitel 5.2.4

Worksheet hydrogen_10.mws

c ITP Bonn                               1995                                          filename: hydrogen.ms

Autor: Komma                                                                             Datum: 31.12.93

Korrigierte und ergänzte Fassung 10.02.04

Index: H-Orbitale, zugeordnete Legendre-Polynome, verallgem. Laguerre-Polynome

Thema: Darstellung der Wahrscheinlichkeitsdichten des H-Elektrons

Synopsis: 3d-Darstellung und Polarplot der Kugelfunktionen, Plot der Radialfunktion,

                Wahrscheinlichkeitsdichte als Funktion von r und theta (3d- und contour-Plot).

                 Demonstration zur Verwendung von Prozeduren und Funktionen.

H-Orbitals

>    restart:

With(orthopoly) macht das package bekannt, in dem die Legendre-Polynome P(l,x) zu finden sind.

>    with(plots): #with(orthopoly);

>    #P(5,x);#test

>    #diff(P(5,x),x$2); # test

>   

Man kann P(l,x) aber auch leicht selbst definieren:

>    P:=(l,x)->if l<>0 then 1/(2^l*l!)*diff((x^2-1)^l,x$l)

>                             else   1 fi;

P := (l, x) -> if l <> 0 then 1/(2^l)/l!*diff((x^2-1)^l,`$`(x,l)) else 1 end if

>    P(5,x); #test # expand(%);

x^5+5*(x^2-1)*x^3+15/8*(x^2-1)^2*x

>    expand(%);

63/8*x^5-35/4*x^3+15/8*x

>   

Dann schreibt man eine kleine Prozedur für die zugeordneten Legendre-Polynome:

>    Y:=proc(l,m,x) ;

>    if m=0 then P(l,x)  else

>    (1-x^2)^(m/2)*diff(P(l,x),x$m)

>    fi;

>    end;

Y := proc (l, m, x) if m = 0 then P(l,x) else (1-x^2)^(1/2*m)*diff(P(l,x),`$`(x,m)) end if end proc

>   

Test

>    Y(5,0,u);

u^5+5*(u^2-1)*u^3+15/8*(u^2-1)^2*u

>    expand(%);

63/8*u^5-35/4*u^3+15/8*u

>   

Für die Plots sollte man x durch cos(theta) ersetzen. Diese Konstruktion ist einfacher, als von vornherein mit theta zu arbeiten (wegen der Ableitungen). Man benötigt aber eine dummy-Variable (ebenfalls wegen diff).

>    ply:=dummy->subs(dummy=cos(theta),Y(l,m,dummy));

ply := dummy -> subs(dummy = cos(theta),Y(l,m,dummy))

>    l:=2:m:=0:

>    ply(poo);

3/2*cos(theta)^2-1/2

Normierung

>    l:='l' : m:='m':

>    Ny2:=(2*l+1)*(l-m)!/(l+m)!; #Normierung^2*4*Pi

Ny2 := (2*l+1)*(l-m)!/(l+m)!

Y(l,m,x);

Kontrolle der "spherical harmonics" nach Vorgabe von l und m:

>    l:=3:m:=2:Y(l,m,x);ply(o);

15*(1-x^2)*x

15*(1-cos(theta)^2)*cos(theta)

>   

>    n:='n';l:='l';m:='m';

n := 'n'

l := 'l'

m := 'm'

>    Ny2*ply(o)^2;

(2*l+1)*(l-m)!/(l+m)!*((1-cos(theta)^2)^(1/2*m))^2/(2^l)^2/l!^2*diff((cos(theta)^2-1)^l,`$`(cos(theta),l+m))^2

>   

>   

Nun kann man sich das Quadrat der  Kugelfunktion zeigen lassen:

>    l:=3:m:=1:
sphereplot(Ny2*ply(o)^2,phi=Pi/2..2*Pi,theta=0..Pi,axes=boxed,scaling=constrained,grid=[15,100]);

[Maple Plot]

>   

Das ist schon ein ganz brauchbares Maß für die Winkelabhängigkeit der Aufenthaltswahrscheinlichkeit eines Elektrons.

Nicht ganz so schön, aber informativ:

>    polarplot(Ny2*ply(o)^2,theta=0..2*Pi,scaling=constrained);

[Maple Plot]

>   

um Pi/2 gedreht zum Vergleich mit der 3d-Darstellung

(und zur Demonstration eines parametrischen Aufrufes):

>    polarplot([Ny2*ply(o)^2,theta+Pi/2,theta=0..2*Pi],scaling=constrained);

[Maple Plot]

>    l:='l': m:='m': n:='n':

Nun die Radialfunktion:

Man benötigt die verallgemeinerten Laguerre-Polynome die im orthopoly-package stehen.

Wer dieses package nicht hat, schreibt:

>    #u=2*r/(a*n)

>    a:='a': # (Bohrscher Radius =1, wird ggf. weiter unten benoetigt)

>    l:='l': m:='m': n:='n':
L:=(j,k,u)->if j<>0 then 1/j!*exp(u)/u^k*diff(u^(j+k)*exp(-u),u$j)

>                               else 1 fi;

L := (j, k, u) -> if j <> 0 then 1/j!*exp(u)/(u^k)*diff(u^(j+k)*exp(-u),`$`(u,j)) else 1 end if

>   

Radialfunktion(u):

>    Ru:=(n,l,u)->u^l*exp(-u/2)*L(n-l-1,2*l+1,u);

Ru := (n, l, u) -> u^l*exp(-1/2*u)*L(n-l-1,2*l+1,u)

Test (l <= n-1)  ... bitte genau hinschauen: l = el, 1 = eins ... oder andere Namen oder Schrifttypen wählen

>    n:=3:l:=1:

>    Ru(n,l,u);

1/u^2*exp(-1/2*u)*exp(u)*(4*u^3*exp(-u)-u^4*exp(-u))

Radialfunktion(r) und Normierung

>    n:='n':l:='l':r:='r':

>    plr:=u->subs(u=2*r/(n*a),simplify(Ru(n,l,u)));

plr := u -> subs(u = 2*r/n/a,simplify(Ru(n,l,u)))

>    Nr2:=4*(n-l-1)!/(n+l)!/(a^3*n^4); #Normierung^2

Nr2 := 4*(n-l-1)!/(n+l)!/a^3/n^4

>    n:=3:l:=1:

>    plr(etwas);

-2/3*exp(-1/3*r/a)*r/a*(-4+2/3*r/a)

>   

>   

Für die Wahrscheinlichkeitsdichte muß die Radialfunktion plr = R(n,l,r) noch mit r multipliziert

werden, dann kann man sich wieder das Quadrat zeigen lassen:

>    a:=1:

>    plot(Nr2*(r*plr(o))^2,r=0..30);

[Maple Plot]

Weitere Tests

>    l:=0:

>    Nr2*(r*plr(o))^2;

1/243*r^2*exp(-1/3*r)^2*(6-4*r+4/9*r^2)^2

>    n:='n';Nr2*(r*plr(o))^2;

n := 'n'

1/(n-1)!/n!/n^2*exp(-r/n)^2*Sum(binomial(n-1,_k1)*pochhammer(n-_k1+1,_k1)*(2*r/n)^(n-_k1)*(-1)^(n-_k1),_k1 = 0 .. n-1)^2

>    plot([seq(Nr2*(r*plr(o))^2,n=1..3)],r=0..30,color=[red,green,blue]);

[Maple Plot]

Test der Normierung

>    seq(int(Nr2*(r*plr(o))^2,r=0..infinity),n=1..3);

1, 1, 1

>    l:=1:

>    plot([seq(Nr2*(r*plr(o))^2,n=2..3)],r=0..30,color=[red,green,blue]);

[Maple Plot]

>    seq(int(Nr2*(r*plr(o))^2,r=0..infinity),n=2..3);

1, 1

>   

>   

>   

Will man noch die Winkelabhängikeit sehen, so stellt man das Produkt (r*R*Y)^2 am besten

dreidimensional und in Polarkoordinaten dar (mit der richtigen nlm-Kombination).

>    n:='n': l:='l':m:='m':

>    plr(o);ply(o);

-(2*r/n)^(-l-1)*exp(-r/n)*Sum(binomial(n-l-1,_k1)*pochhammer(n+l-_k1+1,_k1)*(2*r/n)^(n+l-_k1)*(-1)^(n-l-_k1),_k1 = 0 .. n-l-1)/(n-l-1)!

(1-cos(theta)^2)^(1/2*m)/(2^l)/l!*diff((cos(theta)^2-1)^l,`$`(cos(theta),l+m))

>   

>    Nr2*Ny2*(r*plr(etwas)*ply(o))^2;

4/(n-l-1)!/(n+l)!/n^4*(2*l+1)*(l-m)!/(l+m)!*r^2*((2*r/n)^(-l-1))^2*exp(-r/n)^2*Sum(binomial(n-l-1,_k1)*pochhammer(n+l-_k1+1,_k1)*(2*r/n)^(n+l-_k1)*(-1)^(n-l-_k1),_k1 = 0 .. n-l-1)^2*((1-cos(theta)^2)^(...
4/(n-l-1)!/(n+l)!/n^4*(2*l+1)*(l-m)!/(l+m)!*r^2*((2*r/n)^(-l-1))^2*exp(-r/n)^2*Sum(binomial(n-l-1,_k1)*pochhammer(n+l-_k1+1,_k1)*(2*r/n)^(n+l-_k1)*(-1)^(n-l-_k1),_k1 = 0 .. n-l-1)^2*((1-cos(theta)^2)^(...

>   

>   

>    n:=3:l:=2:m:=1:

>    Nr2*Ny2*(r*plr(etwas)*ply(o))^2;

4/6561*r^6*exp(-1/3*r)^2*(1-cos(theta)^2)*cos(theta)^2

>    plot3d([r*cos(theta),r*sin(theta),Nr2*Ny2*(r*plr(etwas)*ply(o))^2],

>                 r=0..30,theta=Pi/2..2*Pi,axes=boxed);

[Maple Plot]

>   

Nach der Ausführung des obigen Befehls kann man mit der linken Maustaste die box angeklickt halten und so ziehen, daß man das Bild  von oben betrachtet. Der Plot wird nach Betätigen der rechten Maustaste neu aufgebaut. Man wähle dann style=contour!

Oder automatisiert:

>    plot3d([r*cos(theta),r*sin(theta),Nr2*Ny2*(r*plr(o)*ply(o))^2],

>                 r=0..20,theta=0..2*Pi,axes=boxed,orientation=[0,0],

>                 shading=z,style=patchcontour,scaling=constrained);

[Maple Plot]

>   

Mit Titel und ohne Achsen

>    txt:=`nlm=`||n||l||m:

>    txt;

`nlm=321`

>   

>    plot3d([r*cos(theta),r*sin(theta),Nr2*Ny2*(r*plr(o)*ply(o))^2],        r=0..20,theta=0..2*Pi,orientation=[0,0],style=patchcontour,scaling=constrained,shading=z,title=`txt`);

[Maple Plot]

>   

Und weil das so schön ist, macht man eine procedur daraus:

>    myhydrogenplots:=proc(n,l,m) global a,p1,p2,p3,p4,p5; local txt;

>    a:=1;

>    #print(Y(l,m,x));

>    #print(ply(o));

>    txt:=`nlm=`||n||l||m:

>    p1:=sphereplot(Ny2*ply(o)^2,phi=Pi/2..2*Pi,theta=0..Pi,axes=boxed,scaling=constrained,grid=[15,100],title=`txt`);

>    print(p1);

>    p2:=polarplot([Ny2*ply(o)^2,theta+Pi/2,theta=0..2*Pi],scaling=constrained,title=`txt`);

>    print(p2);

>    p3:=plot(Nr2*(r*plr(o))^2,r=0..30,title=`txt`);

>    print(p3);

>    p4:=plot3d([r*cos(theta),r*sin(theta),Nr2*Ny2*(r*plr(o)*ply(o))^2],

>              r=0..30,theta=Pi/2..2*Pi,axes=boxed,title=`txt`);

>    print(p4);

>    p5:=plot3d([r*cos(theta),r*sin(theta),Nr2*Ny2*(r*plr(o)*ply(o))^2],

>                 r=0..30,theta=0..2*Pi,axes=boxed,orientation=[0,0],style=patchcontour,scaling=constrained,

>                  shading=z,title=`txt`);

>    end;

myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...
myhydrogenplots := proc (n, l, m) local txt; global a, p1, p2, p3, p4, p5; a := 1; txt := `nlm=` || n || l || m; p1 := sphereplot(Ny2*ply(o)^2,phi = 1/2*Pi .. 2*Pi,theta = 0 .. Pi,axes = boxed,scaling ...

>   

Bitte warten, bis alle fünf Plots erschienen sind.

>    myhydrogenplots(4,3,0);

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

>   

In der Prozedur wurden die Plots unter Namen abgespeichert und können nun weiter verwendet werden, wenn man z.B. die Optionen ändern will.

>    display(p1,style=wireframe,shading=z,orientation=[56,70]);

[Maple Plot]

>   

Es lohnt sich aber auch das Abspeichern der Plots in .m-files, man kann sich so eine schnell abrufbare Sammlung von "Orbitals" zulegen.

>   

>    for n from 4 to 4 do

>    for l from 0 to n-1 do

>    for m from 0 to l do

>    txt:=`nlm=`||n||l||m:

>    #print(polarplot([Ny2*ply(o)^2,theta+Pi/2,theta=0..2*Pi],scaling=constrained,title=`txt`));

>    p||n||l||m:=plot3d([r*cos(theta),r*sin(theta),Nr2*Ny2*(r*plr(o)*ply(o))^2],

>                 r=0..50,theta=0..2*Pi,axes=boxed,orientation=[0,0],style=contour,scaling=constrained,

>                  shading=z,title=`txt`,numpoints=1000);

>    print(p||n||l||m);

>    od; od; od;

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

>   

Abspeichern

>    #save p4||(0..3)||(0..3),`hydro4.m`;

>   

>    #restart:

>    #read `hydro4.m`;

>    #anames();

>   

>    #p411;

>   

Das Ergebnis kann von jeder Maple-session aus mit

read `hydro4.m`;

eingelesen werden. Mit anames(); erfährt man die Namen der Plots (allerdings auch alle anderen zugewiesenen Namen der jeweiligen session -- anames(); also möglichst zu Beginn der session verwenden).

>   

komma@oe.uni-tuebingen.de