Gross-Pitaevskii equation with NDSolve

Question

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 01 '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 01 '18 at 20:39

Vitaliy Kaurov · Answer 1 · 2018-05-01T21:06:44.273

7

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.

edited May 01 '18 at 21:06

answered May 01 '18 at 20:58

Vitaliy Kaurov

73,078
9
204
355

@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 01 '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 01 '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 01 '18 at 21:39

score 4 · Answer 2 · answered Jan 18 '23 at 08:29

You can run this with FEM, but it will take a bit of time:

tEnd = 1;
Monitor[AbsoluteTiming[
  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, 
      DirichletCondition[u[t, x] == 0, x == -a || x == a]}, 
     u, {t, 0, tEnd}, {x, -a, a}, 
     EvaluationMonitor :> (currentTime = 
        Row[{"t = ", CForm[t], " / ", tEnd}])];], currentTime]

You'd probably also want to use a finer mesh:

Plot3D[Evaluate[Abs[First[u /. sol][t, x]]^2], {t, 0, 1}, {x, -a/1.5, 
  0}, PlotPoints -> 50, ColorFunction -> "Rainbow", PlotRange -> All]

Gross-Pitaevskii equation with NDSolve

2 Answers2

Linked