Note: for messing around with this code without doing a full run, I recommend putting
nframes = 30
gridsize=100
momentum={90,30}
gaussianwidth=0.1
or so. Crank-Nicolson is absolutely stable, but you still only get good results if the wavelength is decently greater than the lattice size.
nframes = 30
gridsize=100
momentum={90,30}
gaussianwidth=0.1
or so. Crank-Nicolson is absolutely stable, but you still only get good results if the wavelength is decently greater than the lattice size.
In[]:=
gridsize=250;L=1.0;dx=L/gridsize;linearAlgebraTime=0;plottingTime=0;videoTime=0;exportTime=0;gaussianwidth=0.1;momentum={230,180};nframes=600;fps=60;dt=0.0001;myModulo[x_]=Mod[x-1,gridsize]+1;index[i_,j_]=gridsize(myModulo[j]-1)+myModulo[i];(*Dirichletconditions*)zeroOrOne[i_,j_]=If[i<1||i>gridsize||j<1||j>gridsize,0.0,1.0];matrixElements=Flatten[Table[{{index[i,j],index[i,j+1]}zeroOrOne[i,j+1],{index[i,j],index[i,j-1]}zeroOrOne[i,j-1],{index[i,j],index[i+1,j]}zeroOrOne[i+1,j],{index[i,j],index[i-1,j]}zeroOrOne[i-1,j],{index[i,j],index[i,j]}-4.0},{i,1,gridsize},{j,1,gridsize}]];hamiltonian=-0.5/dx^2SparseArray[matrixElements,{gridsize^2,gridsize^2}];psi[x_,y_]:=Exp[(-(x-0.5)^2-(y-0.5)^2)/(2gaussianwidth^2)-Imomentum.{x,y}];psivec=Flatten[Table[psi[idx,jdx],{i,1,gridsize},{j,1,gridsize}]];plotPsi:=Module[{t,frame},{t,frame}=AbsoluteTiming[Magnify[Image[Partition[Hue[Arg[#]/(2Pi),1,Clip[2Abs[#]^2,{0,1}]]&/@psivec,gridsize]],3]];plottingTime+=t;frame];(*CrankNicolson.Theimplicitequationis:(psi[i,j,n+1]-psi[i,j,n])/(dt)0.5(hamiltonian.psi[i,j,n+1]+hamiltonian.psi[i,j,n])ThiscanbeformedintoanequationLHS.psi[i,j,n+1]RHS.psi[i,j,n],whichcanbesolvedwithLinearSolvetogetthenexttimestep.*)LHS=IdentityMatrix[gridsize^2,SparseArray]+I0.5hamiltoniandt;RHS=(IdentityMatrix[gridsize^2,SparseArray]-I0.5hamiltoniandt);linsolve=LinearSolve[LHS];step:=Module[{t},t=First@AbsoluteTiming[(psivec=linsolve[RHS.psivec])];linearAlgebraTime+=t;];frames=Table[step;plotPsi,{i,1,nframes}];{videoTime,video}=AbsoluteTiming[Parallelize@FrameListVideo[Rasterize/@frames,RasterSize->300,FrameRate->60]];exportTime=First@AbsoluteTiming[Export["quantumParticle2.mp4",video]];video
Out[]=
In[]:=
Dataset@<|"Linear algebra time (s)"->linearAlgebraTime,"Plotting time (s)"->plottingTime,"mp4 encoding and export (s)"->videoTime+exportTime|>
Out[]=