Temat: Zależne od czasu równanie Schrödingera. Ruch pakietów falowych. Wartości własne.
Przykłady programów w JavaTM .


Na podstawie: Internet

1 Ruch paczki falowej

Poniższy rysunek przedstawia stary sposób ilustrowania zjawiska ruchu paczki falowej w zadanym potencjale (w tym przypadku nad studnią prostokątna).
pic//Sqwell.gif
Ruch pakietu falowego.
Teraz możemy zrobić to lepiej.
Oto, jak np. robi to John L. Richardson z Silicon Graphics.
Oscylacje pakietu falowego w 1dim a tak Denis Rapaport:
1dim
2dim Które przedstawienie jest lepsze? Wydaje się, że ostatnia, chociaż jesli chcemy obejrzeć niektóre szczegóły ...

2 Leapfrog method

Jak możemy to zrealizować sami? Ponieważ nie jesteśmy pierwsi, zobaczmy co zrobili inni. Można skorzystać z metody opracowanej przez Visschera, [P.B. Visscher, Comp. in Phys. , 5 (1991) 596-598], którą można znaleźć także u Horbatscha [Horbatsch, QM]. [The time stepping algorithm is described also in: Richardson, John L., Visualizing quantum scattering on the CM-2 supercomputer, Computer Physics Communications 63 (1991) pp 84-94]. Diabli wiedzą kto był pierwszy. Dobrze jest zajrzeć do klasyków mechaniki kwantowej: L. Schiff, Mechanika kwantowa, wyd. PWN, Warszawa, 1977, str.102 oraz rozdział 12. Schiff cytuje też pracę z 1967r: 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 (1967) 177. Jeszcze jedną interesującą publikacją jest: Galbraith, Ching, Abraham, Am. J. Phys. 52, (1984) 60-68. 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 obliczenia, na zmianę, raz wartości ψR, a raz ψI w kolejnych chwilach t0+∆t/2, t0+∆t), ... w odpowiedni sposób. Z równania Schrödingera mamy mianowicie
ψR′ = H ψI , ψI′ = − H ψR .
Po dyskretyzacji czasu (∆t) i przestrzeni (∆x) odkładamy na osi czasowej odcinki ∆t/2. Ponieważ ewolucja ψR wyznaczana jest przez ψI, natomiast ewolucja ψI wyznaczana jest poprzez ψR to możemy zastosować następujący schemat iteracyjny. Dla części rzeczywistej zapiszemy
ψR(x, t + ∆t) = ψR(x, t) + ∆t H ψI(x, t + ∆t/2)
(1)
a dla urojonej
ψI(x, t + (3/2) ∆t) = ψI(x, t + ∆t/2) − ∆t H ψR(x, t + ∆t) .
(2)
Centralne różnice skończone wyznaczaja pochodne z dokładnością do drugiego rzędu w t. Metoda jest stabilna o ile tylko ∆t < ∆x2. Problem polega tylko na tym, że nie znamy wartości ψR i ψI w tej samej chwili czasu t 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 chwili t, a część urojona ψI powinna być zadana w chwili t+∆t/2.
p(x,t; x0,p0,w0)
=
2
√
w0
 π1/4

√
1 − 2 i t w02
 exp(− p02
(2 w0)2
 )
×exp(w02 (i (x0 − x) − p0/(2 w02))2
(1 − 2 i t w02)
 )
(3)
Schemat programu wygląda jakoś tak. 
/*
 * Scheme of the Wave Packet program
 * (basis: Horbatsch, Maple progarm)
 *
*/
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++) 
    xmesh[i] = -xmax + (i-1)*dx;
// 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;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];
    }
    // Time:=Time+dt;
}
/*
 * 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;
  }
}
.

3 Metoda NCGF rozwiązywania równania Schrödingera

Metoda NCGF (Numerova-Cowella-Goodwina-Foxa) nadaje się do rozwiązywania specjalnych zagadnień własnych. Metoda Numerowa-Cowella
Z a d a n i e 1. Napisz program NCGF w JavaTM . Zastosuj go do rozwiązania równań rozpatrywanych w zadaniach 1-4 zamieszczonych w części Metoda NCGF ....
Jeśli rozpatrujemy rozpad stanów rezonansowych, to przybliżone położenie energii stanów rezonansowych (przy założeniu, że Γ << E) można oszacować rozwiązując równanie Schrödingera metodą NCGF.