Numerické řešení 1-D Schrödingerovy rovnice

Po své první zkoušce z kvantové mechaniky jsem si chtěl ověřit, zda něco z nabytých znalostí dokážu použít k řešení nějakého problému. Vybral jsem si to úplně nejjednodušší: jednorozměrný problém jedné částice. Vše jsem vytvořil bez konzultace s někým, kdo tomu rozumí, nebo na to má alespoň nějaký názor. Následujícímu textu proto nesmíte věřit, berte ho spíše jako inspiraci. Připomínky pište do diskuse.

Řešení jsem hledal pomocí programu Octave. Použil jsem následující diskretizaci

x = linspace (xmin, xmax, nx) .

Pracuji tedy na intervalu (xmin,xmax) rozděleném na nx dílků. Pro vytvoření Hamiltoniánu je nutné vytvořit operátor druhé derivace. Ten je na diskrétní 1-dimenzionální mříži reprezentován následující maticí.

/ -2  1  0  0 ..... 0  1  \
|  1 -2  1  0 ..... 0  0  |
|  0  1 -2  1          .  |
|  .       .           .  |    1
|  .         .         .  |  ------  ,
|  .           .       .  |   dx^2
|  .          1 -2  1  0  |
|  0  0       0  1 -2  1  |
\  1  0 ..... 0  0  1 -2  /
	

kde dx = (xmax-xmin)/nx. Umístil jsem jedničky do rohů, což je zřejmě v souladu periodickými okrajovými podmínkami. Tato matice se v Octavu vytvoří příkazy

diff2matrix = diag(ones(nx-1,1),1) + diag(ones(nx-1,1),-1) + diag(-2*ones(nx,1));
diff2matrix(1,nx)=1;
diff2matrix(nx,1)=1;
diff2matrix /= dx^2;

Dále definuji potenciál, ve kterém se zkoumaná částice bude pohybovat. Volím potenciál jako kvadratickou funkci a řeším tedy lineární harmonický oscilátor.

potencial = x.^2;

Nyní již mohu přistoupit k vytvoření Hamiltoniánu

Ham = -hbar^2/(2*m)*diff2matrix + diag(potencial);

Volím takové jednotky, aby Planckova konstanta i hmotnost částice byly rovny jedné :-). Nalezení vlastních vektorů a vlastních čísel Hamiltoniánu a tedy vyřešení nečasové Schrödingerovy rovnice je nyní již hračkou s užitím funkce eig();

[eigenvec,eigenval] = eig(Ham);

Následující obrázek ukazuje pět vlastních vektorů s nejmenšími vlastními čísly, čili nejnižší energetické stavy. K vlastním funkcím jsem přičetl odpovídající vlastní čísla a zakreslil jsem je do jednoho grafu s potenciálem. Toto běžně používané zobrazení dobře ilustruje vztah vlnových funkcí a potenciálu. Je vidět, že částice se může vyskytnout i v oblasti pod grafem, která je klasickou fyzikou zapovězená.

Obrázek 1.: Vlastní stavy LHO

Abych mohl studovat vývoj obecné vlnové funkce psi_init, musím provést její rozvoj do báze vlastních funkcí Hamiltoniánu. Koeficienty určím následovně:

koef(1:nx)=0;
for i=1:nx
	koef(i) = psi_init*eigenvec(:,i);
endfor

Mohu počítat časový vývoj vlnové funkce

function y = vlna(t)
	y(1:nx)=0;
	for i=1:nx;
		y+=koef(i)*eigenvec(:,i)'*exp(-eigenval(i,i)*t/(I*hbar));
	endfor
endfunction?

A zobrazit ho jako animaci na obrazovce. Viz Obrázek 2.:

for t=linspace(tmin,tmax,nt)
	plot (x,abs(vlna(t)).^2);
endfor

Obrázek 2.: Vlnové klubko v LHO

Zajímavým poznatkem kvantové mechaniky je existence tunelového jevu. Tak se nazývá průchod částice skrz potenciálovou bariéru, vyšší než je energie částice. Namodeloval jsem ukázku tohoto jevu v potenciálu s dvěma prohlubněmi.

Obrázek 3.: Tunelový jev

Kompletní skript v Octavu si můžete stáhnout zde.