Numerical Homework 2

(*Defining descretizing constants*)
Nx=100; xmin=-10.; xmax=10.;dx=((xmax-xmin)/(Nx+1));Nt=500;dt:=((2(π/Ω))/Nt);

(*Exact wavefunction for harmonic oscillator*)
psi[n_,x_]:=(1/(Sqrt[(Sqrt[π])]))(1/(Sqrt[(2^n)n!]))HermiteH[n,x]Exp[-(1/2)(x^2)]//N;

(*Numerical operators: 1D second derivative + position operator*)
(Subscript[d, 2])=(1/(dx^2))SparseArray[{{i_,i_}->-2.,{i_,j_}/;Abs[i-j]== 1->1.},{Nx,Nx}];X=SparseArray[{i_,i_}->xmin+dx i,{Nx,Nx}];

(*General form of 1D Hamiltonian*)
H[t_]:=((-1)/2)(Subscript[d, 2])+V[t]

(*Cayley’s form of time-evolution operator*)
Uplus[t_]:=IdentityMatrix[Nx]+(1/2)I H[t+(dt/2)] dt//N;
Uminus[t_]:=IdentityMatrix[Nx]-(1/2)I H[t+(dt/2)] dt//N;

(*Potential function defintion, with a driving frequency*)
V[t_]:=(1/2)X.X-X*f[t];

(* Simulation protocol*)
Simulate[ω_]:=(Ω=ω;
f[t_]:=Sin[Ω t];
xc[t_]=Simplify[Integrate[f[tp]Sin[t-tp],{tp,0,t}]];Ψ[n_,x_,t_]=Simplify[psi[n,x-xc[t]]Exp[I(-(n+0.5)t+xc'[t](x-((xc[t])/2))+(1/2)Integrate[f[tp]xc[tp],{tp,0,t}])]]//N;Psi=((Eigenvectors[H[0]][[Nx]])/(Sqrt[dx]))//N;data=Table[Psi=LinearSolve[Uplus[(i-1)dt],Uminus[(i-1)dt].Psi];If[Mod[i,(Nt/100)]==0,exactψ[x_]=Ψ[0,x,i dt];exactψ0[x_]=psi[0,x-f[i dt]];Show[{Plot[{Re[exactψ[x]],Im[exactψ[x]],Abs[exactψ[x]],Abs[exactψ0[x]]},{x,-5,5},PlotRange->{-1,1},PlotStyle->{Red,Blue,Black,{Dashed,Green}},AxesLabel->{"x","|Ψ|"},PlotLegends->{(TraditionalForm@Re[ψ]//ToString)<>"_Exact",(TraditionalForm@Im[ψ]//ToString)<>"_Exact",(TraditionalForm@Abs[ψ]//ToString)<>"_Exact",(TraditionalForm@Abs[(Subscript[ψ, 0])]//ToString)<>"_Exact"}],ListPlot[{Re[Psi],Im[Psi],Abs[Psi]},PlotRange->{-1,1},PlotLegends->{(TraditionalForm@Re[ψ]//ToString)<>"_Numerical",(TraditionalForm@Im[ψ]//ToString)<>"_Numerical",(TraditionalForm@Abs[ψ]//ToString)<>"_Numerical"},DataRange->{xmin+dx,xmax-dx},PlotStyle->{Red,Blue,Black}]}],Nothing],{i,1,Nt}];ListAnimate[data,20,AnimationRunning->False])

11.38

f[t_]:=0*Sin[Ω t];Ω=0.5;Psi=(((((Eigenvectors[H[t]][[Nx]])/(Sqrt[dx]))+((Eigenvectors[H[t]][[Nx-1]])/(Sqrt[dx]))))/(Sqrt[2]));Data=Table[Psi=LinearSolve[Uplus[(i-1)dt],Uminus[(i-1)dt].Psi];If[Mod[i,(Nt/100)]==0,exactψ[x_]=(((psi[0,x]Exp[((-I i dt)/2)]-psi[1,x]Exp[((-3I i dt)/2)]))/(Sqrt[2]));Show[{Plot[{Re[exactψ[x]],Im[exactψ[x]],Abs[exactψ[x]]},{x,-5,5},PlotRange->{-1,1},PlotStyle->{Red,Blue,Black},AxesLabel->{"x","|Ψ|"},PlotLegends->{(TraditionalForm@Re[ψ]//ToString)<>"_Exact",(TraditionalForm@Im[ψ]//ToString)<>"_Exact",(TraditionalForm@Abs[ψ]//ToString)<>"_Exact"}],ListPlot[{Re[Psi],Im[Psi],Abs[Psi]},PlotRange->{-1,1},PlotLegends->{(TraditionalForm@Re[ψ]//ToString)<>"_Numerical",(TraditionalForm@Im[ψ]//ToString)<>"_Numerical",(TraditionalForm@Abs[ψ]//ToString)<>"_Numerical"},DataRange->{xmin+dx,xmax-dx},PlotStyle->{Red,Blue,Black}]}],Nothing],{i,1,Nt}];ListAnimate[Data,20,AnimationRunning->False]

11.39 a): Ω=(1/5)ω

Simulate[0.2]

11.39 a): Ω=5ω

Simulate[5.]

11.39 a): Ω=(6/5)ω

Simulate[6./5.]