# 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

## 1 Answer

Advice:

• 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]};

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

Plot3D[Abs[sol[t, x]]^2, {t, 0, 1}, {x, -a/1.5, 0},
PlotPoints -> 50, ColorFunction -> "Rainbow"] 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
• @Gluoncito what I meant is using approximate and asymptotic solutions to derive meaningful spatial and temporal sizes for NDSolve from initial parameters of your equation. There are simple techniques for this classic problem in literature. – Vitaliy Kaurov May 1 '18 at 21:25
• @Vitaly Kaurov, the only thing I can acknowledge is that the initial condition was barely compatible with the boundary condition, but that was not a problem for the algorithm. I'm not an expert in this subject what classical literature could you recommend me? – Gluoncito May 1 '18 at 21:39