Nieskończona studnia 
z barierą w środku

Mamy rozwiązania symetryczne i antysymetryczne. Wystarczy rozważyć prawą
część studni: psi =A sin(k(L-x))
Wewnątrz bariery mamy psi = B1 cosh(kappa x) albo psi = B2 sinh(kappa x)
gdzie k^2= 2 m E /hbar^2 oraz kappa^2 =2 m (V-E) /hbar^2

Warunek zszycia dla rozwiązań symetrycznych

warsym=κ((Sinh[κ a])/(Cosh[κ a]))+k((Cos[k(L-a)])/(Sin[k(L-a)]))

Warunek zszycia dla rozwiązań antysymetrycznych

waranty=κ((Cosh[κ a])/(Sinh[κ a]))+k((Cos[k(L-a)])/(Sin[k(L-a)]))

Pozimy energetyczne dla zadanych parametrów

mamy:

k=((Sqrt[2 m En])/hbar);
κ =((Sqrt[2 m (V-En)])/hbar);

Dalej zakładamy

m=1.0;hbar=1.0;L=2.0;a=1;

warsym

waranty

teraz zakładamy jakąś wartość V i wykreślamy warunki zszycia

V=5;
Plot[{warsym,waranty},{En,0,V},
PlotRange -> {-5, 5},
PlotStyle->{{Blue, AbsoluteThickness[2]},
   {Red, AbsoluteThickness[2]}},
 AxesStyle -> Directive[AbsoluteThickness[1], 16]]

Rozwiązanie numeryczne

Mamy V = 3.

Przybliżone rozwiązania: 2.14 and 2.32, szukamy dokładnych rozwiązań w ich pobliżu:

End1=FindRoot[warsym==0,{En,2.14}]

En1=En/.End1

End2=FindRoot[waranty==0,{En,2.3}]

En2=En/.End2

En2-En1

Funkcje falowe dla dwóch najniższych poziomów

liczby falowe

k1=k/.En->En1;
k2=k/.En->En2;
κ1=κ/.En->En1;
κ2=κ/.En->En2;

Znajdujemy współczynnik A i C w funkcji B

A1=B1 ((Cosh[κ1 a])/(Sin[k1(L-a)]));
A2=B2 ((Sinh[κ2 a])/(Sin[k2(L-a)]));

Nieznormalizowana funkcja falowa dla poziomu podstawowego

psi1[x_]:=If[x<-a,(A1/B1) Sin[k1(x+L)],If[x<a, Cosh[κ1 x],(A1/B1) Sin[k1(L-x)]]]

Plot[psi1[x],{x,-L,L},
PlotStyle->{Red,AbsoluteThickness[2]},
AxesStyle->Directive[AbsoluteThickness[1],16]]

NIeznormalizowana funkcja falowa stanu wzbudzonego

psi2[x_]:=If[x<-a,-(A2/B2) Sin[k2(x+L)],If[x<a,Sinh[κ2 x],(A2/B2) Sin[k2(L-x)]]]

Plot[psi2[x],{x,-L,L},
PlotStyle->{Blue,AbsoluteDashing[5], AbsoluteThickness[2]},
AxesStyle->Directive[AbsoluteThickness[1],16]]

Normalizacja:

N1=NIntegrate[Conjugate[psi1[x]]psi1[x],{x,-L,L}]

N2=NIntegrate[Conjugate[psi2[x]]psi2[x],{x,-L,L}]

Znormalizowane funkcje falowe:

Plot[{((1/(Sqrt[N1]))psi1[x]),((1/(Sqrt[N2]))psi2[x])},{x,-L,L},
PlotStyle->{{Red,AbsoluteThickness[2]},
{Blue,AbsoluteDashing[5], AbsoluteThickness[2]}},
AxesStyle->Directive[AbsoluteThickness[1],16]]

Rozkłady prawdopodobieństwa:

Plot[{((1/(Sqrt[N1]))psi1[x])^2,((1/(Sqrt[N2]))psi2[x])^2},{x,-L,L},
PlotStyle->{{Red,AbsoluteThickness[2]},
{Blue,AbsoluteDashing[5], AbsoluteThickness[2]}},
AxesStyle->Directive[AbsoluteThickness[1],16]]

sprawdzamy ortogonalność:

Chop[NIntegrate[Chop[Conjugate[psi2[x]]psi1[x]],{x,-L,L}]]

Zależność czasowa

Tworzymy zależne od czasu kombinacje linowe zlokalizowane w t=0 po prawej albo lewej stronie bariery

Psi1[t_,x_]:=Exp[-I En1 t](1/(Sqrt[N1]))psi1[x]-Exp[-I En2 t](1/(Sqrt[N2]))psi2[x];
Psi2[t_,x_]:=Exp[-I En1 t](1/(Sqrt[N1]))psi1[x]+Exp[-I En2 t](1/(Sqrt[N2]))psi2[x];

Funkce falowe w t=0:

Plot[{Psi1[0,x],Psi2[0,x]},{x,-L,L},
PlotStyle->{{Red,AbsoluteThickness[2]},
{Blue,AbsoluteDashing[5], AbsoluteThickness[2]}},
AxesStyle->Directive[AbsoluteThickness[1],16]]

Ewolucja czasowa rozkładu prawdopodobieństwa

Animate[
Plot[Conjugate[Psi1[t0,x]]Psi1[t0,x],{x,-L,L},
PlotRange->{0,2.5},
PlotStyle->{Red,AbsoluteThickness[2]},
AxesStyle->Directive[AbsoluteThickness[1],16]],{t0,0,π/(En2-En1)},
AnimationRunning->False]

Animate[
Plot[{Conjugate[Psi1[t0,x]]Psi1[t0,x],Conjugate[Psi2[t0,x]]Psi2[t0,x]},{x,-L,L},
PlotRange->{0,2.5},
PlotStyle->{{Red,AbsoluteThickness[2]},{Blue,AbsoluteThickness[2]}},
AxesStyle->Directive[AbsoluteThickness[1],16]],{t0,0,2π/(En2-En1)},
AnimationRunning->False]