# Computer simulation of Schrödinger equation [duplicate]

I am looking for a computer program which simulates the Schrödinger equation (say for a single particle) in two dimensions and for potentials and initial states specified by the user. Typical situations I have in mind are the double slit experiment, harmonic oscillator, Aharonov–Bohm, tunneling effect etc. Unfortunately I was only able to find the one dimensional case on the web.

So are there some 2d programs available?

PS: I don't think my question is a duplicate of this one. I only ask about a single particle, not about thousands. On the other hand I want freedom in my choice of potential. Also I want a numeric solution of a PDE, while the other question is about (as far as I understand, which is not ver far) a Monte Carlo simulation. Anyway these are completely different requirements.

This probably isn't exactly what you're looking for, but if you're looking for the time-independent bound states of a system, the Fourier grid Hamiltonian method may be applicable. Here is an application of it to the following strange-looking potential well: Here are a few low-energy bound states:       And here are some of the high-energy states:      I used a fairly symmetric-looking potential surface, but you can do arbitrary shapes or squiggly lines if you really wanted to.

Here is Mathematica code used to draw the potential and compute the bound states. Beware that in 2D, just using 81 grid points in each dimension means that you have to diagonalize a $6561\times 6561$ matrix, which can take several minutes. The following defines the spatial axes and constructs the kinetic energy matrix:

\$HistoryLength = 0;
hue = Compile[{{z, _Complex}}, {(1.0 Arg[-z] + \[Pi])/(2 \[Pi]),
Exp[1 - Max[Abs[z], 1]], Min[Abs[z], 1]},
CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
vec[A_] := Flatten[Transpose[A, Reverse[Range[ArrayDepth[A]]]]];
nx1 = 8;
np1 = 8;
nn1 = nx1 np1;
L1 = 8;
\[CapitalDelta]x1 = L1/nx1;
\[CapitalDelta]p1 = 2 \[Pi]/\[CapitalDelta]x1;
\[Alpha]1 = \[CapitalDelta]p1/(2 \[CapitalDelta]x1);
pMax1 = \[CapitalDelta]p1 np1/2;
axis1 = (Range[nn1] - 1/2) L1/nn1;
TRow1 = If[EvenQ[nn1],
1/2 Append[
Table[(2.0 pMax1^2)/nn1^2 (-1)^j/Sin[\[Pi] j/nn1]^2, {j,
nn1 - 1}], pMax1^2/3 (1 + 2.0/nn1^2)],
1/2 Append[
Table[(2.0 pMax1^2)/nn1^2 ((-1)^j Cos[\[Pi] j/nn1])/
Sin[\[Pi] j/nn1]^2, {j, nn1 - 1}],
pMax1^2/3 (1 + 1.0/nn1^2)]];
T1 = MapIndexed[RotateRight, ConstantArray[TRow1, nn1]];
nx2 = 8;
np2 = 8;
nn2 = nx2 np2;
L2 = 8;
\[CapitalDelta]x2 = L2/nx2;
\[CapitalDelta]p2 = 2 \[Pi]/\[CapitalDelta]x2;
\[Alpha]2 = \[CapitalDelta]p2/(2 \[CapitalDelta]x2);
pMax2 = \[CapitalDelta]p2 np2/2;
axis2 = (Range[nn2] - 1/2) L2/nn2;
TRow2 = If[EvenQ[nn2],
1/2 Append[
Table[(2.0 pMax2^2)/nn2^2 (-1)^j/Sin[\[Pi] j/nn2]^2, {j,
nn2 - 1}], pMax2^2/3 (1 + 2.0/nn2^2)],
1/2 Append[
Table[(2.0 pMax2^2)/nn2^2 ((-1)^j Cos[\[Pi] j/nn2])/
Sin[\[Pi] j/nn2]^2, {j, nn2 - 1}],
pMax2^2/3 (1 + 1.0/nn2^2)]];
T2 = MapIndexed[RotateRight, ConstantArray[TRow2, nn2]];


The following allows you to enter a potential surface and plot it (the code below takes Potential to be the strange crosshair-shaped surface I showed at the beginning, but it can be changed to other shapes):

\[Delta] = 0.3;
Potential[x1_,
x2_] := -600 \[Delta]/(Sqrt[(x1 - L1/2)^2] + \[Delta]) \[Delta]/(
Sqrt[(x2 - L2/2)^2] + \[Delta]) -
300 \[Delta]/(
Abs[Sqrt[(x1 - L1/2)^2 + (x2 - L2/2)^2] - 2] + \[Delta]);
Hamiltonian[x1_, p1_, x2_, p2_] := (p1^2 + p2^2)/2 + Potential[x1, x2];
DensityPlot[Hamiltonian[x1, 0, x2, 0], {x1, 0, L1}, {x2, 0, L2},
PlotRange -> All, PlotPoints -> 90, ColorFunction -> GrayLevel]


The following constructs the Hamiltonian and diagonalizes it:

V = DiagonalMatrix[
vec[Table[
Potential[axis1[[k1]], axis2[[k2]]], {k1, nn1}, {k2, nn2}]]];
T = KroneckerProduct[T2, IdentityMatrix[nn1]] +
KroneckerProduct[IdentityMatrix[nn2], T1];
H = T + V;
{\[Lambda], \[CapitalLambda]} = {#1[[Ordering[#1]]], \
#2[[Ordering[#1]]]} & @@ Eigensystem[H];


And the following plots the bound states:

mat[v_] := Partition[v, nn1]\[Transpose];
interp[m_] :=
8 InverseFourier[
RotateLeft[
8 Dimensions[m]], Floor[Dimensions[m]/2]]];
Manipulate[{\[Lambda][[k]],
If[Abs[\[Lambda][[k - 1]] - \[Lambda][[k]]] <= 0.01 ||
Abs[\[Lambda][[k + 1]] - \[Lambda][[k]]] <= 0.01, "Mixed",
"Pure"],
Image[hue[20 interp[mat[\[CapitalLambda][[k]]]]], ColorSpace -> Hue,
Magnification -> 1]}, {k, 1, 500, 1}]


Unfortunately I don't exactly remember what physical units system the code uses, but I do remember that I wrote it partially based on David Tannor's exposition of the FGH method as outlined in his book Quantum Mechanics: A Time-Dependent Perspective.

• Fantastic plots and thanks for the code +1. I'm sure a few people will have fun with that. Have you thought about putting it on the Wolfram Demonstrations site? Also, I love your user page, but is there any chance we could get to know some tiny part of the background of someone who clearly has considerable grasp of many fields of physics? – WetSavannaAnimal Apr 3 '14 at 13:16
• @DumpsterDoofus although your code works (and I have yet to make some sense of it), I get a warning for the first part: "Compile::nogen: A library could not be generated from the compiled function. >>". – Ruslan Apr 6 '14 at 4:49

There are some simulation tools available online, but whether they are useful to you depends on the details of your requirements.

Check out the list of quantum simulators here.
Or this one.

• Thanks! Actually I am looking for something more flexible, where I can play around with different potentials etc. – Jan Weidner Apr 2 '14 at 11:57
• @JanWeidner A very late reply I know, but I recently decided to try making some cool-looking Schrodinger equation simulations, in 2D and 3D. I ended up writing my own app to do the numerical work, and then using ParaView for the rendering. I've posted some of the results on YouTube, and my code and executables on GitHub. The app does let you enter an an arbitrary function for the potential energy. – Paul G Sep 7 '16 at 3:44