Wavefunction time evolution

Time-dependent Harmonic Oscillator

Time-independent harmonic oscillator

division=200.;
mass=2;

DiffMatrix[m_,dx_,Vx_]:=(DiagonalMatrix[-((ℏ^2)/(2m))((-2)/(dx^2))+Vx]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{Length[Vx]-1}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{Length[Vx]-1}],-1])/.ℏ->1

Vfunc[k_,x_]:=(k/2) (x^2)

Varray[V_,boundary_,dx_]:=Table[Vfunc[V,x],{x,-boundary,boundary,dx}];

harmpotential=Varray[60,5,2 5/division];

Set up the Matrix for the Problem and calculate the Eigenfunctions and Eigenvectors:

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
answers=Eigensystem[MProblem];
{evals,evects}=Sort[Transpose[answers],#1[[1]]<#2[[1]]&]//Transpose;
stateenergies=Take[evals,10]

Plot the Eigenvectors on top of the Eigenenergies :

potplot=ListPlot[harmpotential,Filling->1.1Min[harmpotential],Joined->True,Frame->True,PlotRange->{0,60},GridLines->{None,stateenergies},GridLinesStyle->Directive[Thickness[0.002],Red],Axes->False];
ScaleFactor=10;
stateplot=ListPlot[Table[ScaleFactor evects[[i]]+evals[[i]],{i,1,10}],Joined->True,Filling->Table[n->evals[[n]],{n,1,10}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False];
Show[stateplot,potplot,PlotRange->{Automatic,{0,60}}]

Time-dependent harmonic oscillator :

Frequency :

ω=(Sqrt[D[Vfunc[60,x],{x,2}]/mass]);
tfinal=2π/ω;

Start in the groundstate :

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
uinitial=evects[[1]];

Energy of the ground state :

{uinitial.MProblem.uinitial,N[ω/2]}

(Subscript[Ω, 0])=(ω/2) in the ground state (vacuum energy)

NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}]

answertdepground[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}]

ListPlot[Abs[answertdepground[tfinal]],PlotRange->All]

Realteil der Wellenfunktion :

Animate[Show[potplot,ListPlot[ evals[[1]]+ScaleFactor Re@answertdepground[t],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{-5,30}],{t,0,2tfinal},AnimationRunning->False]

Absolutwert der Wellenfunktion (Aufenthaltswahrscheinlichkeit):

Animate[Show[potplot,ListPlot[ evals[[1]]+ScaleFactor Abs@answertdepground[t],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{-5,30}],{t,0,2tfinal},AnimationRunning->False]

Start in the 1st excited state :

expected frequency (Subscript[Ω, 1])=((3ω)/2)

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
uinitial=evects[[2]];
answertdep1[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}];

Realteil der Wellenfunktion :

Animate[Show[potplot,ListPlot[ evals[[2]]+20 Re[answertdep1[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{0,30}],{t,0,2tfinal},AnimationRunning->False]

Absolutwert der Wellenfunktion (Aufenthaltswahrscheinlichkeit):

Animate[Show[potplot,ListPlot[ evals[[2]]+20 Abs[answertdep1[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{0,30}],{t,0,2tfinal},AnimationRunning->False]

Start in the 2nd excited state :

(Subscript[Ω, 2])=((5ω)/2)

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
uinitial=evects[[3]];
answertdep2[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}];

Realteil der Wellenfunktion :

Animate[Show[potplot,ListPlot[{evals[[3]]+20 Re[answertdep2[t]],evals[[2]]+20 Re[answertdep1[t]],evals[[1]]+20 Re[answertdepground[t]]},PlotStyle->{Black,Red,Green}],PlotRange->{-5,30}],{t,0,2tfinal},AnimationRunning->False]

Absolutwert der Wellenfunktion (Aufenthaltswahrscheinlichkeit):

Animate[Show[potplot,ListPlot[{evals[[3]]+20 Abs[answertdep2[t]],evals[[2]]+20 Abs[answertdep1[t]],evals[[1]]+20 Abs[answertdepground[t]]},PlotStyle->{Black,Red,Green}],PlotRange->{-5,30}],{t,0,2tfinal},AnimationRunning->False]

Start in superposition state :

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
uinitial=1/(Sqrt[2])(evects[[2]]+evects[[1]]);
answertdepsup[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}];

Animate[Show[potplot,ListPlot[ evals[[1]]+ScaleFactor Abs[answertdepsup[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{50,150},{0,20}}],{t,0,2tfinal},AnimationRate->0.4,AnimationRunning->False]

Particles bounces in between the potential walls.

Start in a localized (Gaussian) state:

PDF[NormalDistribution[0,0.05],x]

ugauss=Table[PDF[NormalDistribution[0,0.05],x],{x,-5,5,2 5/division}];

Normalize to |⟨Ψ|Ψ⟩(|^2)=1:

normsqrt=(Sqrt[ugauss\[Conjugate].ugauss])

uinitial=ugauss/normsqrt;

ListPlot[uinitial,PlotRange->All,Joined->True,Frame->True,FrameLabel->{"x","|Ψ(x)((|)^2)"},BaseStyle->{FontSize->14}]

Energy of the gaussian state:

gaussenergy=uinitial\[Conjugate].MProblem.uinitial

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
answertdeploc[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,tfinal}];

Animate[Show[potplot,ListPlot[gaussenergy+10Abs[answertdeploc[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{60,150},{0,60}}],{t,0,tfinal},AnimationRunning->False]

Start in a localized state at x≠ 0 :

uinitial=#/Norm[#]&[Table[PDF[NormalDistribution[-.3,0.4],x],{x,-5,5,2 5/division}]];
ListPlot[uinitial,PlotRange->All]

Energy of the displaced gaussian state:

gaussenergy2=uinitial\[Conjugate].MProblem.uinitial

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
answertdeploc2[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,tfinal}];

Animate[Show[potplot,ListPlot[gaussenergy2+10Abs[answertdeploc2[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{50,150},{4,12}}],{t,0,tfinal},AnimationRunning->False]