# Gross-Pitaevskii equation with NDSolve

I want to solve the Gross-Pitaevskii equation with NDSolve, I tried the following code:

a = 100;
sig = a/10;
sol = NDSolve[{I D[u[t, x], t] == (-1/2) D[u[t, x], {x, 2}] +
x^2/2 u[t, x] - Abs[u[t, x]]^2 u[t, x],
u[0., x] == Exp[-((x + 30)^2./(2*sig^2))], u[t, a] == 0,
u[t, -a] == 0}, u, {t, 0, 1130}, {x, -a, a}, MaxStepSize -> 0.05,
AccuracyGoal -> 4, PrecisionGoal -> 4];

Animate[Plot[Evaluate[Abs[u[t, x] /. First[sol]]^2], {x, -a, a},
PlotRange -> {0, 10}], {t, 0, 413}]

DynamicNDSolve::ndsz: At t == 0.9724326016597982, step size is effectively
zero; singularity or stiff system suspected.

DynamicNDSolve::eerr: Warning: scaled local spatial error estimate of
1420.198482935478 at t = 0.9724326016597982 in the direction of
independent variable x is much greater than the prescribed error tolerance.
Grid spacing with 4001 points may be too large to achieve the desired
accuracy or precision. A singularity may have formed or a smaller grid
spacing can be specified using the MaxStepSize or MinPoints method options.


I also tried setting the methods:

Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral"}}, AccuracyGoal -> 4,
PrecisionGoal -> 4]


without success.

Reformulation: as recommended below I read some papers of the literature, and I tried to reproduce their results with Mathematica. Everything goes fine until I added rotation to the BEC. I tried the following code:

a = 8;
sig = a/10;
Epsilon = 1.;
Kappa = 1;
Omega_s = 0.;
tfin = 2 Pi;
sol = NDSolve[{I Epsilon D[u[t, x, y],
t] == (-Epsilon^2/2) (
D[u[t, x, y], {x, 2}] + D[u[t, x, y], {y, 2}]) + (x^2 + y^2)/
2 u[t, x, y] +
I Epsilon  Omega_s (x D[u[t, x, y], {y, 1}] -
y D[u[t, x, y], {x, 1}]) + Kappa  Abs[u[t, x, y]]^2 u[t,
x, y], u[0.,x,y] == (1/Sqrt[Pi sig] ) Exp[-(x^2 + y^2)/(2 sig)],
u[t, a, y] == u[t, -a, y], u[t, x, a] == u[t, x, -a]},
u, {t, 0, tfin}, {x, -a, a}, {y, -a, a}];

Animate[Plot[Evaluate[Abs[u[t, x, 0] /. First[sol]]^2], {x, -a, a},
PlotRange -> {0, 3}], {t, 0, tfin}]


which runs fine if the rotation frequency Omega_s is close to zero, but when I try to reach the critical value for vortex formation (Omega_s=0.25) the code break.

I already did it using the pseudospectral method described above!, very cool Mathematica ! (I'll try to apply the FEM package to see if it can resolve this problem as well).

• I tried removing the AccuracyGoal, PrecisionGoal, and MaxStepSize, and limiting the integration to $(0,1)$, which should reproduce your error according to the value of $t$ in the error. However, it worked fine, albeit slowly. Do you know that you need such a tiny step size over such a long time period? You could also use MaxStepFraction to ensure that a certain number of steps are taken. – MarcoB May 1 '18 at 19:42
• @MarcoB Yes, it works removing the goals, I didn't expect it. I was able to extend the integration time as well. Thank you. Please note my edit in the second part of the question. – Gluoncito May 1 '18 at 20:39

• Use cyclic (looped) boundary conditions at the end points if GP wave function is the same there: u[t,a]==u[t,-a]

• Try to make sense of your physical parameters not to run the solver into computational death

• Use only exact symbolic integers and rationals when setting up equations inside NDSolve, do not use decimals of limited precision

• NDSolveValue is better for direct plotting

Here are simple ranges where it works (start from here and modify reading tutorials on proper settings):

a=50;
sig=a/10;

eqs={I D[u[t,x],t]==(-1/2) D[u[t,x],{x,2}]+x^2/2 u[t,x]-Abs[u[t,x]]^2 u[t,x],
u[0,x]==Exp[-((x+15)^2/(2 sig^2))],u[t,a]==u[t,-a]};


sol = NDSolveValue[eqs, u, {t, 0, 1}, {x, -a, a}];

Second part sounds as a separate question. Generally use Interpolation or NonlinearModelFit or BSplineFunction` to turn discrete data into a smooth functional surface for IC.
• @Vitaly Kaurov, thanks, I already simulated it with the suggestion of MarcoB. About your comment: "Try to make sense of your physical parameters not to run the solver into computational death", does not $t>1$ make physical sense to you?, or anything else? – Gluoncito May 1 '18 at 21:11