Následující příklad ukazuje, že v Mathematice přímá aplikace elementárních rovnic kvatové mechaniky umožňuje velmi zhruba načrtnout evoluci stavu v 1D potenciálové jámě.

Samozřejmě, jsou známé efektivnější a přesnější metody.

ϕ[i_] :=x^i Exp[-x^2/2]

Takto bázové funkce vypadají:

Table[ ϕ[i],{i,0,3}]
Plot[%,{x,-5,5},PlotRange->Full]

Skalaární součin v 1D je obyčejný integrál, pro jednoduchost bázových funkcí asi vystačíme s analytickýni výpočty

dot[f_]:=Integrate[ f, {x,-Infinity,Infinity}]

Můžeme si dovolit jen konečný počet bázových funkcí, takže půjde o přibližný výpočet

dim = 5; (*časem zvýšíme*)

Připomeneme si pořadím v jakém Table pracuje se dvěma cykly

Table[Subscript[a,{i,j}],{i,0,dim},{j,0,dim}]//MatrixForm

matU = ParallelTable[dot[ϕ[i]*ϕ[j]],{i,0,dim},{j,0,dim}]

opH[f_] := -D[f,{x,2}]+x^2 f

matH = ParallelTable[dot[opH[ ϕ[i] ]*ϕ[j]],{i,0,dim},{j,0,dim}]

Eigenvalues[ Inverse[matU] . matH]

Zdá se, že jsme dostali rozumné hodnoty pro vlastní hodnoty energie harmonického oscilátoru, tak zkusíme jiný potenciál

opH[f_] := -D[f,{x,2}]+x^4 f

matH = ParallelTable[dot[opH[ ϕ[i] ]*ϕ[j]],{i,0,dim},{j,0,dim}];

{dim, Eigenvalues[ Inverse[matU] . matH]//N}

dim = 11;
matU = ParallelTable[dot[ϕ[i]*ϕ[j]],{i,0,dim},{j,0,dim}];
matH = ParallelTable[dot[opH[ ϕ[i] ]*ϕ[j]],{i,0,dim},{j,0,dim}];
{dim, Eigenvalues[ Inverse[matU] . matH]//N}

dim = 14;
matU = ParallelTable[dot[ϕ[i]*ϕ[j]],{i,0,dim},{j,0,dim}];
matH = ParallelTable[dot[opH[ ϕ[i] ]*ϕ[j]],{i,0,dim},{j,0,dim}];
{dim, Eigenvalues[ Inverse[matU] . matH]//N}

Pozorování: nejnižších 5 vlastních hdodnot energie by snad mohlo být dobře

λ = Eigenvalues[{matH,matU}//N]

ev = Eigenvectors[ Inverse[matU] . matH//N]//Chop;

basis = Table[ϕ[j],{j,0,dim}]

Vlastní funkce získáme vynásobením bázových funkcí příslušnými koeficienty

ev[[-1]].basis (*základní stav*)
Plot[- %, {x,-4,4}]

-ev[[-4;;-1]].basis//Simplify
Plot[ %, {x,-4,4},PlotRange->Full,PlotTheme->"Detailed",PlotLegends->Table[i,{i,3,0,-1}]]

Řešení časové Schroedingerovy rovnice

Nyní budeme bykreslovat řešení časové Schroedingerovy rovnice
Vezmeme nějaký počáteční stav, rozložíme ho do vlastních vektorů hamiltoniánu a následně využujeme faktu, že ve vlastní bázi H je evoluční operátor diagonální.

vlastni = ev.basis//Reverse;

psi0 = Exp[-2(x+1)^2]
Plot[ psi0,{x,-4,4}]

a0 = ParallelTable[ dot[psi0 vlastni[[i]]]/dot[ vlastni[[i]]^2],{i,1,dim+1}];

{psi0, a0.vlastni};
Plot[ %,{x,-6,6},PlotRange->Full]

λ  ... vlastni hodnoty,
vlastni ... senznam funkci x pro jednotlive i = 1 .. dim + 1
a0 ... slozky poc . stavu ve vlastni bazi H 

napsieme
Ψ[t] funkce v case t v promenne x

Ψ[t_Real]:= a0 . DiagonalMatrix[Exp[- I t λ ]] . vlastni

{psi0^2, Abs[Ψ[1.0]]^2};
Plot[ %,{x,-6,6},PlotRange->Full]

Animate[ 
Plot[{psi0^2, Abs[Ψ[t]]^2},{x,-6,6},PlotRange->Full],
{t,0,0.1}
]