# How can I solve a three-dimensional Gross-Pitaevskii equation?

I want to solve a three-dimensional Gross-Pitaevskii equation with an uniform potential. The related papers are Emergence of a Turbulent Cascade in a Quantum Gas and Synthetic dissipation and cascade fluxes in a turbulent quantum gas.

The code I tried was

hbar = 1;

kB = 1.380649*10^-23;
\[Mu] = 1;(*2 10^-9 kB;*)
\[Xi] = 1;(*1.2 10^-6;*)

Udepth = 30 \[Mu];
Uz = 27/1.2 \[Xi];
UR = 16/1.2 \[Xi];

N0 = 1. 10^5;
n0 = N0/(Pi UR^2 Uz);
g = \[Mu]/n0;
\[CapitalDelta]U = \[Mu]; (*shaking amplitude*)
\[Omega]shaking =
2 Pi 8;

Lx = 40 \[Xi]; Ly = 40 \[Xi]; Lz = 40 \[Xi];
x0 = Lx;
M = 1;
Vtrap[x_, y_, z_, t_] :=
If[x^2 + y^2 < UR && -Uz/2 < z < Uz/2, 0,
Udepth] + \[CapitalDelta]U Sin[\[Omega]shaking t] z/Uz
eqs = { D[\[Psi]1[t, x, y, z],
t] == ((-g n0 Laplacian[\[Psi]1[t, x, y, z], {x, y,
z}]) + (Vtrap[x, y, z, t] +
g (Abs[\[Psi]1[t, x, y, z]]^2)) \[Psi]1[t, x, y, z])/(
I g n0 hbar)};
bc = {\[Psi]1[0, x, y, z] ==
Exp[-x^2/(10 \[Xi])^2] Exp[-y^2/(10 \[Xi])^2] Exp[-z^2/(10 \
\[Xi])^2](*1/Sqrt[rcx rcy Pi]Exp[-x^2/rcx^2/2]Exp[-y^2/rcy^2/2]*),
DirichletCondition[\[Psi]1[t, x, y, z] == 0, True]};
solv = NDSolve[{eqs, bc}, \[Psi]1, {t, 0, 0.5}, {x, -Lx/2,
Lx/2}, {y, -Ly/2, Ly/2}, {z, -Lz/2, Lz/2},
Method -> {"PDEDiscretization" -> {"MethodOfLines", \
{"SpatialDiscretization" -> "FiniteElement"}}}]


But this costs too much time and can not output. What should I do?

• You expect a real solution \[Psi]1? Mar 6 at 8:58
• @UlrichNeumann No. The wavefunction can be a complex one. Mar 6 at 9:58

After proper normalization we have

hbar = 1;

kB = 1.380649*10^-23;
\[Mu] = 1;(*2 10^-9 kB;*)
\[Xi] = 1;(*1.2 10^-6;*)

Udepth = 30 \[Mu];
Uz = 27/1.2 \[Xi];
UR = 16/1.2 \[Xi];

N0 = 1. 10^5;
n0 = N0/(Pi UR^2 Uz);
g = \[Mu]/n0;
\[CapitalDelta]U = \[Mu]; (*shaking amplitude*)
\[Omega]shaking = 2 Pi 8;

Lx = 40 \[Xi]; Ly = 40 \[Xi]; Lz = 40 \[Xi];
x0 = Lx; g1 = g/(g n0 );
M = 1;
Vtrap[x_, y_, z_, t_] :=
If[x^2 + y^2 < UR && -Uz/2 < z < Uz/2, 0,
Udepth] + \[CapitalDelta]U Sin[\[Omega]shaking t] z/(Uz g n0)
eqs = {I hbar D[\[Psi]1[t, x, y, z],
t] == ((-Laplacian[\[Psi]1[t, x, y, z], {x, y, z}]) + (Vtrap[x,
y, z, t] + g1 (Abs[\[Psi]1[t, x, y, z]]^2)) \[Psi]1[t, x, y,
z])};
bc = {\[Psi]1[0, x, y, z] ==
Exp[-x^2/(10 \[Xi])^2] Exp[-y^2/(10 \[Xi])^2] Exp[-z^2/(10 \
\[Xi])^2](*1/Sqrt[rcx rcy Pi]Exp[-x^2/rcx^2/2]Exp[-y^2/rcy^2/2]*),
DirichletCondition[\[Psi]1[t, x, y, z] == 0, True]};

solv = NDSolve[{eqs, bc}, \[Psi]1, {t, 0, 0.5}, {x, -Lx/2,
Lx/2}, {y, -Ly/2, Ly/2}, {z, -Lz/2, Lz/2}]


Visualization

Table[Plot3D[
Evaluate[Abs[\[Psi]1[t, 0, y, z] /. solv[]]], {y, -Ly/2,
Ly/2}, {z, -Lz/2, Lz/2}, PlotLabel -> t,
ColorFunction -> "ArmyColors", Mesh -> None, PlotRange -> All], {t,
0, .5, .1}] • Thanks a lot! @Alex And can I ask you how much time this code takes in your computer? Mar 6 at 13:02
• Ah, I see that time is about 5771.71 s. Did you ask about time optimization? Mar 6 at 17:49
• Yeah,@Alex. After a whole night, the updated code is still running in my computer(Macbook pro 13.3' M1, 16G+512G ). So what is the possible reason here? Is this issue very dependent on the performance of the computer? And can I ask you how much memory do you use? Mar 7 at 2:57
• And the vesion is mathematica 13.0. Mar 7 at 3:32
• May be solution depends on version? In my case \$Version is 13.2.1 for Microsoft Windows (64-bit) (January 27, 2023). Machine has 32 GiB physical memory, and 33.85 GiB virtual memory. Mar 7 at 3:35