I am trying to follow the main ideas presented in this question, applying it to my own problem, which is a complex, time-dependent, nonlinear PDE:

$$i \frac{\partial \psi}{\partial t} = \left[ -\nabla^2 + (|\psi|^2-1) \right] \psi$$

It is supposed to show the evolution of a vortex, so that is what I take as initial condition. The code is organized as follows

I work in 2D, on a regular grid from xl to xr and yl to yr, on the time interval [0.,tmax].

xl = yl = -5.;
xr = yr = 5.;
tmax = 2.;

This is the phase of the (complex) initial condition:

ϕ[x_?MachineNumberQ, y_?MachineNumberQ] := ArcTan[x, y]
ϕ[0., 0.] := 0.

Taking the gradient of the phase and plotting it with StreamDensityPlot, one can see the winding velocity of the "quantum" fluid.


As for the spatial part of the initial condition, I take a Tanh profile, because the density $|\psi|^2$ in the center of the vortex should be 0, and far away it should be 1 (the equation in normalized in that way).

vortex[x_?MachineNumberQ, y_?MachineNumberQ] := 
Tanh[4 Norm[{x, y}]] E^(I ϕ[x, y])

Here is the plot of the spatial profile: spatial-profile

Now for the main part of the code, which uses NDSolve to propagate the initial condition in time, with the appropriate boundary conditions. Finally, the solution is plotted at different instances of time and an animation is created.

 sol = ψ /. 
           {(-1. + Abs[ψ[x, y, t]]^2)*ψ[x, y, t] - I*Derivative[0, 0, 1][ψ][x, y, t] - 
              Derivative[0, 2, 0][ψ][x, y, t] - Derivative[2, 0, 0][ψ][x, y, t] == 0, 
              ψ[x, y, 0.] == vortex[x, y], 
              ψ[xl, y, t] == 1., 
              ψ[xr, y, t] == 1., 
              ψ[x, yl, t] == 1., 
              ψ[x, yr, t] == 1.
            }, ψ, {x, xl, xr}, {y, yl, yr}, {t, 0., tmax}, 
            Method -> {"MethodOfLines", 
                       "DifferentiateBoundaryConditions" -> False,
                       "SpatialDiscretization" -> 
                             "DifferenceOrder" -> 4, 
                             "MaxPoints" -> 100, 
                             "MinPoints" -> 100, 
                             "AccuracyGoal" -> 2, 
                             "PrecisionGoal" -> 2}

Export["evolve.gif", pl, AnimationRepetitions -> Infinity, "DisplayDurations" -> .1]

There are various errors in the output of this code. Apart from the (I hope) benign complaint about the inconsistency of the boundary and initial condition, NDSolve simply chokes at t=0.02, complaining about convergency issues. I have tried playing with the parameters, but I have not found any way of circumventing this issue.


Here is an example:

showStatus[status_] := 
clearStatus[] := showStatus[""];

xl = yl = -5;
xr = yr = 5;
tmax = 2;

ϕ[x_, y_] := ArcTan[x, y]

vortex[x_, y_] := Tanh[4 Norm[{x, y}]] E^(I ϕ[x, y])

eqn = I*Derivative[0, 0, 1][ψ][x, y, 
     t] == -Laplacian[ψ[x, y, t], {x, 
       y}] + (Abs[ψ[x, y, t]]^2 - 1)*ψ[x, y, t];

(* boundary and initial conditions should not be inconsistent. This message is not benign *)

bcs = {ψ[xl, y, t] == vortex[xl, y], ψ[xr, y, t] == 
    vortex[xr, y], ψ[x, yl, t] == 
    vortex[x, yl], ψ[x, yr, t] == vortex[x, yr]};
ics = ψ[x, y, 0] == vortex[x, y];

(* the key point here is to use the "Pseudospectral" difference order *)
nxy = 33;
sol = NDSolveValue[{eqn, ics, bcs}, ψ, {x, xl, xr}, {y, yl, 
   yr}, {t, 0, tmax}
  , Method -> {"MethodOfLines", "SpatialDiscretization" -> {
      , "MaxPoints" -> nxy, "MinPoints" -> nxy
      , "DifferenceOrder" -> "Pseudospectral"
  , EvaluationMonitor :> showStatus["t = " <> ToString[CForm[t]]]

Animate[Plot3D[Im[sol[x, y, t]], {x, xl, xr}, {y, yl, yr}], {t, 0, 
  tmax}, AnimationRunning -> False]

enter image description here

You can find more information about the "Pseudospectral" difference order in tutorial/NDSolvePDE. Hope this helps.

  • $\begingroup$ @Kuba, thanks for the edit! $\endgroup$
    – user21
    Jul 9 '13 at 10:05
  • 1
    $\begingroup$ Your wellcome :), since You are posting a lot of differential-equations etc., You may want to look at this halirutan's script $\endgroup$
    – Kuba
    Jul 9 '13 at 10:10
  • $\begingroup$ The "pseudospectral" difference order, from what I understand from the docs, applies FFT methods, which basically assumes periodic boundary conditions, right? Whereas mine were fixed (Dirichlet type). $\endgroup$
    – Andrei
    Jul 9 '13 at 12:19
  • $\begingroup$ @Andrei, I think the first thing you should look at is the boundary/initial conditions. Those need to be consistent. Then one could look further how to solve the problem. $\endgroup$
    – user21
    Jul 9 '13 at 12:25
  • 2
    $\begingroup$ "Simulating as if it were infinite" is a very difficult task in general. It may be easiest to just extend the spatial domain so that boundary effects don't kick in until after tmax... $\endgroup$
    – Jens
    Jul 9 '13 at 18:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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