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).

  • 2
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    May 1, 2018 at 19:42
  • $\begingroup$ @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. $\endgroup$
    – Gluoncito
    May 1, 2018 at 20:39

2 Answers 2



  • 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):


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


enter image description here

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"]

enter image description here

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.

  • $\begingroup$ @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? $\endgroup$
    – Gluoncito
    May 1, 2018 at 21:11
  • $\begingroup$ @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. $\endgroup$ May 1, 2018 at 21:25
  • $\begingroup$ @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? $\endgroup$
    – Gluoncito
    May 1, 2018 at 21:39

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

tEnd = 1;
  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]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.