A. Baran
Kwantowa paczka falowa
Spis treści
Spis rysunków . . . 1
Ruch paczki falowej . . . 2
Leapfrog method . . . 2
Metoda Numerova . . . 5
Unitarność . . . 6
Postać Cayley’a operatora ewolucji . . . 6
Metody stabilne, unitarne i jawne . . . 6
Równanie dyfuzji . . . 8
Więcej wymiarów . . . 8
Problemy graficzne . . . 8
Rzutowanie . . . 9
Poziomice . . . 9
Problem niewidocznych linii . . . 10
Model przestrzeni kolorów RGB . . . 11
Rzutowanie z R2 na przestrzeń RGB . . . 11
Paczka falowa 2D . . . 14
Literatura . . . 20
Spis rysunków
1. Ruch pakietu falowego. . . 22. R . . . 9
3. Interpolacjaz w punkcie wewnętrznym (xm; ym)trójkąta (i-j-k). . . 9
4. Prostopadłościan. . . 10
5. Sześcian kolorów RGB i wartości parametrówr; g; b w jego wierzchołkach . . . 11
6. Sześcian kolorów RGB i wartości parametrów r, g, b w jego wierzchołkach . . . 11
7. Sześcian kolorów RGB. . . 12
8. Rzutowanie R2 (C) na S2 . . . 12
9. Dwa stożki zamiast sfery S2 . . . 13
Ruch paczki falowej
Zależne od czasu równanie Schrödingera. Ruch pakietów falowych. War-tości własne.Przykłady programów w JavaT M.
Poniższy rysunek przedstawia stary sposób ilustrowania zjawiska ruchu paczki falowej w zadanym potencjale (w tym przypadku nad studnią prostokątna).
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 8 9 10 11 START Probability P(x,t) x
Wave packet motion through the Square-well potential.
Rysunek 1. Ruch pakietu falowego.
Obecnie, gdy komputery są szybkie i posiadają olbrzymie możliwości gra-ficzne, możemy zrobić to dużo lepiej.
Tak robi to John L. Richardson z Silicon Graphics. http://www.sgi.com/ fun/java/john/wave-sim.html. Oscylacje pakietu falowego – 1-dim.
A tak Denis Rapaport: http://www.ph.biu.ac.il/~rapaport/java-apps/ quant1d.html. – 1dim; http://www.ph.biu.ac.il/~rapaport/java-apps/ quant2d.html – 2-dim.
Które przedstawienie jest lepsze? Wydaje się, że Rapaporta.
Leapfrog method
Jak możemy to zrealizować? Ponieważ nie jesteśmy pierwsi, zobaczmy co <zrobili inni.
Można skorzystać z metody opracowanej przez Visschera [3], którą można znaleźć także u Horbatscha [4, 6]. Nie do końca wiadomo kto był pierwszy. Dobrze jest zajrzeć do klasyków mechaniki kwantowej [1]. Schiff cytuje też pracę Goldberga i Scheya z 1967r: Jeszcze jedną interesującą publikacją jest praca [5].
Metoda nosi nazwę leap frog czyli ıżabiego skoku. Polega ona na tym, że zespolona funkcja falowa = R+ I, jest zadana w chwili t0 i jej ewolucja, zgodnie z czasowym równaniem Schrödingera, wyznaczana jest poprzez obli-czenia, na zmianę, raz wartości R, a raz I w kolejnych chwilacht0+t=2, t0+t), ... w odpowiedni sposób.
Z równania Schrödingera mamy mianowicie
0
R=H I; I0 = −H R:
Po dyskretyzacji czasu (t) i przestrzeni (x) odkładamy na osi czasowej od-cinkit=2. Ponieważ ewolucja Rwyznaczana jest przez I, natomiast
ewolu-cja Iwyznaczana jest poprzez Rto możemy zastosować następujący schemat iteracyjny. Dla części rzeczywistej zapiszemy
R(x; t + t) = R(x; t) + tH I(x; t + t=2) (1)
a dla urojonej
I(x; t + (3=2)t) = I(x; t + t=2) − tH R(x; t + t) : (2) Centralne różnice skończone wyznaczaja pochodne z dokładnością do drugiego rzędu wt. Metoda jest stabilna o ile tylko t < 2x.
Problem polega tylko na tym, że nie znamy wartości R i I w tej samej chwili czasut lecz naprzemiennie w co drugim kroku o długości t=2. Można je jednak wyznaczyć w tej samej chwili przez interpolację. Dotyczy to również warunku początkowego. Rzeczywista część funkcji falowej R jest zadana w chwilit, a część urojona I powinna być zadana w chwilit + t=2.
p(x; t; x0; p0; w0) = 2 pw 01=4 p 1 − 2itw2 0 exp(− p 2 0 (2w0)2) exp(w2 0 (i(x0−x) − p0=(2w20))2 (1 − 2itw20) ) (3) Schemat programu pokazano poniżej.
/*
* Scheme of the Wave Packet program * (basis: Horbatsch, Maple program) *
*/
static double dt, dth, Time, dxsq; static int nx;
static double[] xmesh, psiR, psiI, Vpot, Hpsi;
// time discretization // time step
dt = 1d/50d;
// half time step dth = dt/2d; // space discretization // space step dx = 1d/5d; dxsq = 1d/(2*dx^2); // space extension xmax = 15d;
// number of space points nx = (int)2*xmax/dx + 1;
// x mesh
for(int i=0;i<nx;i++)
// initialize real and imaginary parts of Psi // at t = 0 and t = dt/2 respectively w0 = 1d/2d; x0 = -10d; p0 = 2d; for(int i=0;i<nx;i++) {
psiR[i] = pack(xmesh[i], 0,x0,p0,w0).Re(); // at t= 0 psiI[i] = pack(xmesh[i],dth,x0,p0,w0).Im(); // at t=dth }
// picture psiR[i], psiI[i]
// and the density |psiR[i]+psiI[i]|^2 // ...
// Hamiltonian on the mesh
// Kinetic energy couples nearest neighbours // Potential energy is local
// Potential; square well for(int i=0;i<nx;i++)
Vpot[i] = Math.abs(xmesh[i]) <= 2 ? -4 : 0;
// picture Vpot[i]
// Calculate H psiR(0), H psiI(dth) doHpsi(psiR);
doHpsi(PsiI);
// doHpsi
public void doHpsi(double[] args) { // calculates H psiR or HpsiI for(int i=1;i<nx-1;i++) Hpsi[i] = -dxsq * (args[1][i+1]+args[1][i-1]-2*args[1][i]) +Vpot[i]*args[1][i]; Hpsi[0] = -dxisq*(args[1][2]-2*args[1][1])+Vpot[1]*args[1][1]; Hpsi[nx] = -dxisq*(args[1][nx-1]-2*args[1][nx]) +Vpot[nx]*args[1][nx]; } /*
* The basic time evolution step for the real and imaginary * parts of the wavefunction. It assumes that the imaginary * part is a half-step ahead of the real part.
*/
public void timeStep() { for(int i=0;i<nx-1;i++) {
psiR[i] = psiR[i] + dt*Hpsi[i]; psiI[i] = psiI[i] - dt*Hpsi[i]; }
}
/*
* Set up the time-loop. */
while ( go ) { timeStep();
// plot xmesh[i],psiR[i]^2+psiI[i]^2],psiR[i],psiI }
// "go" is true if the process is running and false otherwise
class complex { double re, im;
complex() { re = 0; im = 0; }
complex(double x, double y){ re = x;
im = y; }
public void add(complex a, complex b){ re = a.re + b.re;
im = a.im + b.im; }
public void mult(complex a, complex b){ re = a.re * b.re - a.im * b.im; im = a.re * b.im + a.im * b.re; }
public void set(complex a){ re = a.re; im = a.im; } } .
Metoda Numerova
Metoda Numerova-Cowella-Goodwina-Foxa nadaje się do rozwiązywania spe-cjalnych zagadnień własnych, równań, gdzie nie występuje pierwsza pochodna.
Z a d a n i e 1.
Zakoduj algorytm Numerova w JavaT M. Zastosuj program do rozwiązania równań rozpatrywanych w zadaniach 1-4 zamieszczonych w części Metoda Numerova . . . .
Jeśli rozpatrujemy rozpad stanów rezonansowych, to przybliżone położenie energii stanów rezonansowych, przy założeniu, że szerokość rezonansu << E, można oszacować rozwiązując równanie Schrödingera metodą Numerowa.
Unitarność
Zależne od czasu równanie Schrödingera. Ruch pakietów falowych. Sche-mat jawny, stabilny i unitarny. [?].
Pokazany w poprzednim wykładzie schemat iteracyjny rozwiązania równa-nia Schrödingera zależnego od czasu jest schematem jawnym (tzn. n+1 jest wyliczane na podstawie poprzedniego kroku n) lecz nie jest schematem unitar-nym. Operator ewolucji exp(it ^H) zastąpiliśmy tam wyrażeniem 1 it ^H, które nie zapewnia unitarności funkcji falowej (R| |2dx 6= 1). Podobnie jest z prostymi schematami niejawnymi opartymi na tym wyrażeniu. Eliminują one niestabilność numeryczną lecz problem niezachowania unitarności wciąż pozostaje. Przyjrzyjmy się stabilnej metodzie niejawnej. Otrzymamy ją na-stępująco. Funkcja falową n+1w chwilit potraktujmy chwilowo jako daną. Na jej podstawie wyliczmy n. Ponieważ (t) = exp(−itH) (0)}, więc (0) = exp(it ^H) (t). Stąd, rozwijając w szereg Taylora, mamy
n
j = n+1j − (it=x2)[ j+1n+1−2 jn+1− j−1n+1−x2Vj n+1j ]:
Równanie to, chociaż stabilne, daje niejawnie n+1. Metoda jest, jak to nazy-wamy, typu implicit. Aby wyznaczyć n+1potrzeba stosować specjalne techniki algebraiczne.
Ważniejszy jest jednak problem unitarności . Postać Cayley’a operatora ewolucji
Prostym przybliżeniem operatora ewolucji exp(it ^H) jest postać podana przez Cayleya
(1 −1
2t ^H)=(1 + 1
2t ^H) ; (4)
która dodatkowo zapewnia dokładność do rzędut2.
Stosując postać Cayleya operatora ewolucji możemy teraz zapisać
n+1
j = (1 − 12t ^H)=(1 +12t ^H) jn; (5)
Po prostych obliczeniach dostaniemy następujacy schemat różnicowy
n+1 j+1 + (i2x 2 t −t2Vj−2) jn+1+ n+1j−1 = − nj+1+ (i2x 2 t +t2Vj+2) nj − j−1n : (6) Równanie to jest stabilne, zapewnia unitarność, lecz nie daje jawnie wartości dla n+1.
Metody stabilne, unitarne i jawne
Opieramy się na pracy [6].Dyskutowane poprzednio metody były oparte na przybliżeniu dla operatora ewolucji funkcji falowej, który jest formalnym rozwiązaniem równania Schrödin-gera. Równanie ewolucji jest postaci
(x; t) = exp(−it − t0 ¯h
^
H) (x; t0): (7) Najprostsza metoda (niestabilna) była oparta na przybliżeniu
exp(−it ¯h ^ H) = 1 − it ¯h ^ H + O([t ¯h ^ H]2): (8)
Prowadziło to do rozbieżności metody.
Drugie przybliżenie było oparte na postaci Cayleya operatora ewolucji (patrz 4). Ta metoda była wprawdzie stabilna numerycznie lecz niejawna i kosztowna numerycznie gdyż wymagała złożonych obliczeń prowadzących do rozwikłania niejawnego schematu ewolucji.
Metoda, którą podaje Richardson, jest również oparta na przybliżeniu ope-ratora ewolucji, lecz jest to metoda zarówno jawna, jak i unitarna, i stabilna numerycznie.
Operator ewolucji zapisywany jest najpierw w postaci sumy operatora ener-gii kinetycznej ^T i potencjalnej ^V . Dalej ^T rozbija się na sumę M częsci
^
T = PMl=1T^l, dla których łatwo jest policzyć exp ^Tl. Następnie korzystamy z tzw. przybliżenia iloczynowego Trottera:
exp(−it ¯h ^ H) = 0 @YM j=1 exp(−it ¯ h H)^ 1 A exp(−it¯h V ) + O(t^ 2) (9)
Błąd jaki tutaj popełniamy zależy od wartości komutatorów [Tl; Tl0], [Tl; V ] i od wielkości kroku czasowegot.
W przypadku jednego wymiaru przestrzennego podział operatora ^T można określić z formuły różnicowej dla Laplasjanu:
( ^T )n= ¯h
2
2ma2[2 n− n−1− n+1]:
MacierzT jest postaci
[T ] = ¯h 2 2ma2 0 B B B B B B B @ 2 −1 0 : : : 0 −1 −1 2 −1 : : : 0 0 −1 2 −1 : : : 0 ::: 0 : : : −1 2 −1 −1 0 : : : 0 −1 2 1 C C C C C C C A (10)
Macierz ta daje się zapisać jako suma macierzy
[Te] = ¯h 2 2ma2 0 B B B B @ M 0 0 : : : 0 0 M 0 : : : 0 ::: 0 : : : 0 M 0 0 0 : : : 0 M 1 C C C C A (11) oraz [To] = ¯h 2 2ma2 0 B B B B B B B @ 1 0 0 : : : −1 0 M 0 : : : 0 0 0 M : : : 0 ::: 0 : : : 0 M 0 −1 0 : : : 0 1 1 C C C C C C C A ; (12)
gdzieM jest macierzą kwadratową 2 2 M = 1 −1 −1 1 : (13)
Macierze [Te] i [To] są sumami prostymi macierzy M = ¯h2
2ma2M. Łatwo jest obliczyć exponent macierzyM. Mamy mianowicie
exp(−it ¯ h M) = 1 + [exp(−i) − 1] ma2 ¯ h2 M = 1 2 1 + exp(−i) 1 − exp(−i) 1 − exp(−i) 1 + exp(−i) = (14)
gdzie = t¯h=ma2, = [1 + exp(−i)]=2, = [1 − exp(−i)]=2. Macierz exp(−it¯hM) jest unitarna i posiada wartości własne 1=1, 2=exp(−i). Z a d a n i e 2.
Proszę sprawdzić słuszność ostatniej formuły. Obliczyć wartości własneM i exp(−iM).
Zapiszmy wyniki w postaci operatorowej
( ^Te )n= ¯h 2 2ma2 n−12(1 − (−1)n) n−1−12(1 + (−1)n) n+1 ; (15) ( ^To )n= ¯h 2 2ma2 n−12(1 + (−1)n) n−1−12(1 − (−1)n) n+1 : (16) Exponent tych wyrażeń jest
exp(−it=¯h ^Te ) n= n− 2(1 − (−1)n) n−1− 2(1 + (−1)n) n+1 ; (17) exp(−it=¯h ^To ) n= n− 2(1 + (−1)n) n−1− 2(1 − (−1)n) n+1 : (18) Równanie dyfuzji
Wyniki dotyczące stabilności, dyskutowane tutaj dla przypadku równania Schrödingera, są również słuszne dla równania dyfuzji. Podstawieniet → −it iV = 0 przekształca wszystkie wzory otrzymane dla równania Schrödingera w klasyczne równanie dyfuzji. Operatory unitarne zamieniane sa rzeczywistymi operatorami o wartościach własnych mniejszych niż jeden i stabilność metody jest zapewniona.
Więcej wymiarów
W przypadku większej liczby wymiarów wyrażenia są niewiele bardziej skom-plikowane.
Z a d a n i e 3.
Otrzymać wzory różnicowe dyskutowanej metody dla przypadkun-wymiarowego równania Schrödin-gera lub równania dyfuzji.
Problemy graficzne
Ponieważ dyskutowane wyżej problemy można przedstawić graficznie zaj-miemy się przez chwilę sposobami reprezentacji danych na ekranie. Dyskusja będzie dotyczyć trzech spraw. Pierwsza to proste rzutowanie figur 3 wymiaro-wych, druga to rysowanie poziomic (konturów), a trzecia to problem niewidocz-nych linii na rysowaniewidocz-nych rzutach.
X Y Z C P(X,Y,Z) P' p(x,y) obiekt rzut o O x y oko E RZUTOWANIE P'' P'' E O p'' o E P' y x y=Y(Z - Z )/(Z - Z)E o E x=X(Z - Z )/(Z - Z)E o E p' p'' p' Rysunek 2. R Rzutowanie
Łatwo sprawdzić słuszność następujących równości
x = |OP0| |EO|=|EQ| = X(Z
E−ZO)=(ZE−Z) (19)
y = Y (ZE−ZO)=(ZE−Z) (20)
Z a d a n i e 4.
Napisz program w JavaT M(klasa), który rzutuje figurę o zadanych wierzchołkach (xi; yi; zi)na ekran. Zastosuj do kilku wybranych figur przestrzennych.
Poziomice
Poziomice są liniami jednakowej wartości funkcji z = f(x; y), a więc okre-ślone są równaniemz = const. Najwygodniej jest podzielić obszar określoności funkcjif na trójkątne podobszary. Interpolacja wielkości z w takim elementar-nym obszarze trójkąta przebiega następująco (patrz rysunek).
i j k (xi,yi,zi) (xm,yu,zu) (xm,ym,zm) (xm, yl, zl)
Rysunek 3. Interpolacjaz w punkcie wewnętrznym (xm; ym)trójkąta (i-j-k).
Najpierw znajdujemy najdalej wysunięte punkty wx: xmin=min(xi; xj; xk); xmax=max(xi; xj; xk):
W celu zbadania całego obszaru trójkąta przesuwamy linię x = xm od xmin doxmax. Dla każdej wartościxmznajdujemy granice zmiennościy. Nazwiemy
je yl oraz yu. Są one ograniczone dwoma bokami trójkąta i − j − k. Trzeci bok nie jest ważny i wykluczymy go z rozważań. Aby sprawdzić, który to bok sprawdzamy czy xm leży poza przedziałem (xp; xq), gdzie p; q = i; j; k. Jeśli np. xmleży poza przedziałem (xk; xj)(patrz rysunek), to bokij wykluczamy z procedury znajdowaniayu i yl dlax = xm. Jeślixm przekroczyxk to znajdzie się na zewnątrz (xi; xk)i w konsekwencji wykluczymy bok ik. Jeśli będziemy wyznaczaćxu ixl z równań bokówik oraz ij, jak na rysunku, to dostaniemy
yu =yi+ (yk−yi)(xm−xi)=(xk−xi); yl=yi+ (yj−yi)(xm−xi)=(xj−xi):
Jeżelixk =xi lubxj =xi to yu iyl są równe maksymalnej i minimalnej war-tościy. Podobne formuły interpolacyjne stosujemy dla z otrzymując zu izlw punktach (xm; yu)i (xm; yl)(wystarczy zamienicy-ki przez z-ty w równanich). Znajomość graniczl i zu na linii xm pozwala znaleźć wartość z w punkcie ym przez interpolację
z(xm; ym) =zm=zl+ (ym−yl)(zu−zl)=(yu−yl):
Zmieniającxm w granicach (xmin; xmax iymoraz ym w granicach (yl; yu) mo-żemy znaleźćz we wszystkich interesujących nas punktach trójkąta i − j − k. Należy przy tym dobrać wielkość przyrostówx, y tak by otrzymać w miarę gładkie krzywe.
Z a d a n i e 5.
(Dla zainteresowanych) Napisz program dokonujący triangulacji dziedzinyz = f(x; y) i znajdujący poziomicez = const.
Problem niewidocznych linii
Które linie rysowanej na ekranie figury są widoczne, a które nie? Jak je odróżnić? Problem ten prosto opisał H.N. Lerman [7]. Aby zrozumieć na czym to polega spójrzmy na rysunek prostopadłościanu.
1 2 3 4 5 6 7 8 Rysunek 4. Prostopadłościan.
Zewnętrzne wierzchołki prostopadłościanu wraz z łączącymi je krawędziami są zawsze widoczne. Wewnętrzne wierzchołki mogą być albo widoczne albo niewidoczne. Zależy to od tego, ź której strony kartki papieruóne się znajduje, a więc zależy to od wartości współrzędnejz wierzchołka. Np., jeśli kartkę papieru (ekran) ustawimy tak by przechodziła przez z = 0 to wierzchołki o z < 0 (i punkty linii oz < 0 ) nie powinny pojawić się na rysunku.
Z a d a n i e 6.
(Dla zainteresowanych) Napisz program, który rysuje powierzchnięz = f(x; y) bez linii niewidzial-nych.
Model przestrzeni kolorów RGB
Kolorowanie i przedstawianie zespolonej funkcji falowej w przestrzeni kolorów RGB [8, 6].
Model barw RGB (Red, Green, Blue) jest najczęściej stosowanym standar-dem przestrzeni barw. Jest to przestrzeń zwarta i taka, że punkt (rgb)=(000) odpowiada w niej czerni, a punkt (111) bieli. Na przekątnej (000)-(111) sze-ścianu RGB leżą odcienie szarości. Innym wierzchołkom szesze-ścianu odpowiadaja barwy podstawowe: czerwona (010), zielona (100) oraz niebieska (001). Po-zostałe barwy wierzchołkowe to żółta (110) oraz magenta (011). Pośrednie wartości natężeń rgb dają w wyniku mieszankę barw, a więc barwę złożoną.
red green yellow black blue cyan white magenta
Rysunek 5. Sześcian kolorów RGB i wartości parametrówr; g; b w jego wierzchołkach
Standard RGB stosowany jest w urządzeniach wyświetlających takich, jak monitory oraz lampy kineskopowe TV. Obok modelu RGB stosuje się często mo-del barw HSB (Hue, Saturation, Brightnes). Inne standardy, takie jak CMYK (Cyan, Magenta, Yellow) stosowane są w drukarstwie .
Rysunek 6. Sześcian kolorów RGB i wartości parametrów r, g, b w jego wierzchołkach
Dalej pokażemy, jak odwzorować w przestrzeń RGB funkcję określoną na zespolonej płaszczyźnie C (lub R2) tak by można było śledzić jej wartości oraz fazy. Metoda ta nadaje się np. do przedstawienia rozpraszania kwantowej paczki falowej na dwuwymiarowym potencjale (patrz Richardson).
Rzutowanie z R2 na przestrzeń RGB
W celu przedstawienia liczby zespolonej (x + iy) lub pola wektorowego z R2 zastosujemy pewne odwzorowanie C w przestrzeń [0; 1]3. Jedną z wielu możliwości jest odwzorowanie stereograficzne płaszczyzny na sferę jednostkową S2, polegające na tym, że promienie przechodzą przez biegun północny (0; 0; 1)
(000) (100) (011) (010) (110) (001) (111) (101) RGB CUBE r g b
Rysunek 7. Sześcian kolorów RGB
sfery oraz punkty (x; y; 0) na płaszczyżnie C. Punkty z wnętrza koła jednost-kowego odwzorowywane są na półkulę południową, punkt z = 0 na biegun południowy, a punkty leżące poza okręgiem jednostkowym są odwzorowane na półkulę północną. Wszystkim punktom w nieskończoności odpowiada punkt bieguna północnego (0; 0; 1). Okrąg jednostkowy przechodzi w tym odwzoro-waniu w równik sfery. Sytuację przedstawia schematycznie rysunek.
(x, y) (0,0,1)
(0,0,-1) (x',y',z')
Rysunek 8. Rzutowanie R2 (C) na S2
Opisane odwzorowanie jest postaci (patrz również rysunek)
(x0; y0; z0) = 2x; 2y; x
2+y2−1
x2+y2+1 : (21)
Następny krok polega na zanurzeniu sfery jednostkowej w przestrzeń kolo-rów. W tym celu zmniejszymy promień sfery do1=2 i umieścimy jej środek w punkcie (r; g; b) = (1=2; 1=2; 1=2). Oś południe-północ sfery utożsamimy z kie-runkiem przekątnej sześcianu RGB (000) − (111), a więc z szarym kierunkiem czerń-biel w przestrzeni barw. W ten sposób biegun południowy, bliski punk-towi (0; 0; 0) w RGB, jest prawie czerny - brak amplitudy, a biegun północny prawie biały - maksymalna amplituda. Aby określić jednoznacznie orientację azymutalną zażądamy by kierunek osiy przebiegał przez maksimum natężenia czerwieni, a więc przez punkt (010) w przestrzeni RGB.
Opisane odwzorowanie dobrze oddaje sytuację dla punktów z sąsiedztwa okręgu jednostkowego. Początek układu oraz punkty w nieskończoności nie są dostatecznie różne. Aby poprawić sytuację, można zamienić sferę parą stożków o wierzchołkach w punktach (0; 0; 0) i (1; 1; 1), sklejonych na równiku sfery. W ten sposób nie zepsujemy za bardzo obszaru równikowego i jednocześnie "do-damy więcej punktów"w obszarze biegunów sfery, a więc łatwiej będzie rozróżnić punkty w okolicy początku układu i w nieskończoności (rysunek).
(0,0,0)
(1,1,1)
Rysunek 9. Dwa stożki zamiast sfery S2
Odwzorowanie to ma postać
(x0; y0; z0) = 2x; 2y; sign(z)(x
2+y2−1)p3(x2+y2+1 − 2px2+y2)
x2+y2+1 (22)
Dzieląc przez 2 i obracając tak by osie stożków pokryły się z kierunkiem diagonali sześcianu kolorów oraz przesuwając do środka sześcianu kolorów otrzy-mamy ostatecznie r = 12 +sign(x2+y2−1)(1 2 − p x2+y2 x2+y2+1) + 2x p 6(1 + x2+y2) (23) g = 12+sign(x2+y2−1)(1 2− p x2+y2 x2+y2+1) − x p 6(1 + x2+y2)+ y p 6(1 + x2+y2) (24) b =12+sign(x2+y2−1)(1 2− p x2+y2 x2+y2+1) − x p 6(1 + x2+y2)− y p 6(1 + x2+y2) (25) Przykład tego typu odwzorowania pokazuje rysunek. Przedstawiona funkcja f(x; y) została tylko pomalowana, tzn. (x; y) → (r; g; b) z użyciem algorytmu Richardsona bez specjalnych powodów.
Program w JavaT M, który wyznacza kolory podany jest poniżej. \import java.awt.Color;
...
private Color toRGB(float x, float y){
// transform (x,y,z) position to color space rgb float xy, sq, iq, eta, q1, q2, q3, r, g, b; final float half = 0.5f;
Color color; xy = x*x+y*y;
sq = (float)Math.sqrt(xy); iq = 1f/(xy+1f);
Rysunek 10. Funkcja pokolorowana wg. algorytmu Richardsona q1 = half+eta*(half-sq*iq); q2 = iq/(float)Math.sqrt(6); q3 = iq/(float)Math.sqrt(2); r = q1 + 2*x*q2; g = q1 - x*q2 + y*q3; b = q1 - x*q2 - y*q3; color = new Color(r, g, b); return color;
}
Paczka falowa 2D
/*** Scattering of the two dimensional wave packet * * @author AB, 2001 * * @see JL Richardson, CPC 63, 1991, 84-94. * */ import Jampack.*; import Jampack.Times; import java.awt.Color;
public class WavePacket {
static int xpts[], ypts[], nx, ny, drx, dry; final static Znum CZERO = new Znum(0, 0);
double tmin, t, dt;
double x0, xmin, xmax, dx, xscale; double y0, ymin, ymax, dy, yscale;
double hbar, mass, epsilon, xwidth, vx, ywidth, vy; double vwidth, energy, energyScale;
Znum alpha, beta; Zmat psi, expv;
Zmat betaeven, betaodd; double [][] re, im;
WavePacket() { xmin = -4.; xmax = 4.; ymin = -4.; ymax = 4.; nx = 200; ny = nx; xpts = new int[nx]; ypts = new int[ny]; dx = (xmax-xmin)/(nx-1); dy = (ymax-ymin)/(ny-1); xscale = (nx-0.5)/(xmax-xmin); yscale = (ny-0.5)/(ymax-ymin); drx = (int)(xscale*dx + 1); dry = (int)(yscale*dy + 1); re = new double[nx][ny]; im = new double[nx][ny]; }
public void go(){ // set physics energyScale = 1; try {
Physics();
} catch (JampackException e){
System.out.println("Some errors...Jampack"); } StdDraw.enableDoubleBuffering(); StdDraw.setCanvasSize(2*nx,2*ny); StdDraw.setXscale(0,2*nx); StdDraw.setYscale(0,2*ny); while(true){ t += dt/10; // .15 StdDraw.clear(StdDraw.BLACK); doStep(2); makeGraph(); StdDraw.show(); StdDraw.pause(50); //System.out.println(t+" "+dt); //if(Math.round(t*100) % 60==0) //StdDraw.save("pac-"+Math.round(t) +".jpg"); } } /**
* Physical data of the wave packet */
private void Physics() throws JampackException { // This is like in J.L. Richardson, CPC 63, 84-94, 1990
x0 = -1.2; y0 = -1.2;
tmin = 0; t = tmin; hbar = 1; mass = 100; xwidth = 0.4; ywidth = 0.4; vwidth = 0.6; vx = .6; vy = .6; double ds2 = dx*dx+dy*dy; dt = 0.5 * mass * ds2 / hbar; epsilon = hbar * dt /( mass * ds2); double sin = 0.5*Math.sin(epsilon); double cos = 0.5*Math.cos(epsilon); alpha = new Znum(0.5 + cos,-sin); beta = new Znum(0.5 - cos,-sin);
energy = 0.5 * mass * (vx * vx + vy * vy);
// initialize coordinates and the packet double ex, xval, ey, yval, ang, r; for(int x=0;x<nx;x++){
xval = xmin + dx*x;
xpts[x] = (int)(xscale*(xval - xmin));
ex = Math.exp(-((xval-x0)/xwidth)*((xval-x0)/xwidth)); for(int y=0;y<ny;y++){
yval = ymin + dy*y;
ypts[y] = (int)(yscale*(yval - ymin));
ey = Math.exp(-((yval-y0)/ywidth)*((yval-y0)/ywidth)); r = ex*ey; ang = mass*(vx*xval+vy*yval)/hbar; re[x][y] = r*Math.cos(ang); im[x][y] = r*Math.sin(ang); } }
psi = new Zmat(re, im); psi = Star.o(norm(psi), psi);
// init expv matrix for(int x=0;x<nx;x++){
xval = xmin + dx*x; for(int y=0;y<ny;y++){
yval = ymin + dy*y;
r = Potential(xval, yval)*dt/hbar; re[x][y] = Math.cos(r);
im[x][y] = -Math.sin(r); }
}
expv = new Zmat(re, im);
// betaeven(), betaodd() are defined to be zero on odd // or even sites respectively and beta otherwise betaeven = new Zmat(nx, ny);
betaodd = new Zmat(nx, ny); for(int x=0;x<nx;x++)
for(int y=0;y<ny;y++){ if((x+y)%2!=0){ betaeven.put0(x, y, CZERO); betaodd.put0(x, y, beta); } else { betaeven.put0(x, y, beta); betaodd.put0(x, y, CZERO); } } } /** * Potential function * @param x coordinate * @param y coordinate * @return the Potential */
double Potential(double x, double y){ //wall or row
//return x+y < 0 ? 0 : energy*energyScale; // jama kolowa lub pien
return Math.sqrt(x*x + y*y) < vwidth ? -energy*energyScale : 0; //return chpot(x, y);
//return 0; }
double chpot(double x, double y) { // 1/cosh()
double ar = Math.sqrt((x*x + y*y))/vwidth; double ch = 2./(Math.exp(-ar)+Math.exp(ar)); ch = ch*ch;
return -.01*ch*energy*energyScale; }
/**
* The norm of the wavefunction */
private double norm(Zmat psi){ Znum c = new Znum();
double d = 0; for(int x=0;x<nx;x++) for(int y=0;y<ny;y++){ c = psi.get0(x, y); d += c.re*c.re + c.im*c.im; } return 1/Math.sqrt(d*dx*dy); } /**
* Stepping method for two dimensional * time dependent Schroedinger equation */
private void doStep(int num) { try {
for(int k=0;k<num;k++){ psi = Star.o(expv, psi);
psi = Plus.o(Star.o(alpha, psi),
Plus.o(Star.o(betaeven, Cshift.o(psi, 1, -1)), Star.o(betaodd, Cshift.o(psi, 1, +1)))); psi = Plus.o(Star.o(alpha, psi),
Plus.o(Star.o(betaeven, Cshift.o(psi, 1, +1)), Star.o(betaodd, Cshift.o(psi, 1, -1)))); psi = Plus.o(Star.o(alpha, psi),
Plus.o(Star.o(betaeven, Cshift.o(psi, 2, -1)), Star.o(betaodd, Cshift.o(psi, 2, +1)))); psi = Plus.o(Star.o(alpha, psi),
Plus.o(Star.o(betaeven, Cshift.o(psi, 2, +1)), Star.o(betaodd, Cshift.o(psi, 2, -1)))); // norm
//System.out.println("norm = "+norm(psi)); psi = Star.o(norm(psi), psi);
}
} catch (JampackException e) {} }
/**
Makes graph of <b>dens, all, re, im</b></br> depending on no = 1, 2, 3, 4 respectively */
private void makeGraph(){ double f, g;
double amp = 5; Znum p = new Znum(); for(int n=1;n<5;n++){ for(int x=0;x<nx;x++) for(int y=0;y<ny;y++){ p = psi.get0(x, y); switch (n) { case 1: // density f = (p.re*p.re+p.im*p.im)*amp; StdDraw.setPenColor(toRGB(f, f/2)); StdDraw.filledRectangle(xpts[x], ny+ypts[y], drx/2, dry/2); break;
case 2: // real part f = p.re*amp; g = (p.re*p.re+p.im*p.im)*amp; StdDraw.setPenColor(toRGB(f, g)); StdDraw.filledRectangle(nx+xpts[x], ny+ypts[y], drx/2, dry/2); break;
case 3: // imaginary part f = p.im*amp; g = (p.re*p.re+p.im*p.im)*amp; StdDraw.setPenColor(toRGB(g, f)); StdDraw.filledRectangle(xpts[x], ypts[y], drx/2, dry/2); break;
case 4: // both re and im StdDraw.setPenColor(toRGB(p.re*amp, p.im*amp)); StdDraw.filledRectangle(nx+xpts[x], ypts[y], drx/2, dry/2); break; } } } StdDraw.setPenColor(StdDraw.WHITE); StdDraw.text(30,2*ny-20,"abs(f)"); StdDraw.text(nx+30,2*ny-20,"Re(f)"); StdDraw.text(30,nx-20,"Im(f)"); StdDraw.setPenColor(StdDraw.BOOK_BLUE); StdDraw.line(0, ny, 2*nx, ny);
StdDraw.line(nx, 0, nx, 2*ny); } /** * Calculates RGB color; (x, y) -> to (r, g, b) * @param x * @param y - coordinates on C2 * @return Color */
private static Color toRGB(double x, double y){ // transforms (x, y) to color space rgb double xy, sq, iq, eta, q1, q2, q3, r, g, b; final double half = 0.5;
Color color; xy = x*x+y*y; sq = Math.sqrt(xy); iq = 1/(xy+1); eta = xy-1<=0?-1:1; q1 = half+eta*(half-sq*iq); q2 = iq/Math.sqrt(6); q3 = iq/Math.sqrt(2); r = q1 + 2*x*q2; g = q1 - x*q2 + y*q3; b = q1 - x*q2 - y*q3; r = r<0 ? 0 : r; g = g<0 ? 0 : g; b = b<0 ? 0 : b;
color = new Color((float)r, (float)g, (float)b); return color;
}
public static void main(String[] argv){ WavePacket pac = new WavePacket(); pac.go();
} }
Literatura
[1] L. Schiff, Mechanika kwantowa. PWN, Warszawa, 1977, str.102 oraz rozdział 12. [2] A. Goldberg, H.M. Schey. J.L. Schwartz, Computer-Generated Motion Pictures of One-Dimensional Quantum Mechanical Transmission and Reflection Phenomena, Am. J. Phys, 35, 177 (1967).
[3] P.B. Visscher, Comp. in Phys., 5, 596-598 (1991).
[4] M. Horbatsch. Quantum Mechanics Using Maple. Springer, 1995.
[5] Galbraith, I., Ching, Y. S., & Abraham, E. American Journal of Physics, 52, 60 (1984).
[6] Richardson, John L., Visualizing quantum scattering on the CM-2 supercompu-ter, Computer Physics Communications 63, 84-94 (1991).
[7] H.N. Lerman. Software Age, July, 1970.
