Complex valued 2+1D nonlinear PDE using NDSolve

Question

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:

ClearAll[ϕ]
ϕ[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.

currents

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

ClearAll[vortex];
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 = ψ /. 
       First[
         NDSolve[
           {(-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" -> 
                            {"TensorProductGrid", 
                             "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.

score 9 · Accepted Answer · edited Jul 09 '13 at 09:21

9

Here is an example:

showStatus[status_] := 
  LinkWrite[$ParentLink, 
   SetNotebookStatusLine[FrontEnd`EvaluationNotebook[], 
    ToString[status]]];
clearStatus[] := showStatus[""];
clearStatus[]

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

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

ClearAll[vortex];
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" -> {
      "TensorProductGrid"
      , "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.

edited Jul 09 '13 at 09:21

Kuba

136,707
13
279
740

answered Jul 09 '13 at 09:09

@Kuba, thanks for the edit! – Jul 09 '13 at 10:05
1

Your wellcome :), since You are posting a lot of [tag:differential-equations] etc., You may want to look at this halirutan's script – Kuba Jul 09 '13 at 10:10
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). – Andrei Jul 09 '13 at 12:19
@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. – Jul 09 '13 at 12:25
Trying you code, NDSolve chokes with the error "Maximum number of 10000 steps reached at the point t == 1.92". Also, could you please edit your answer to plot Abs[sol]^2 instead of the imaginary part please? That is the physical quantity of interest ;) – Andrei Jul 09 '13 at 12:38
About the boundary conditions, what I can say based on physical arguments is that far away from the center (at infinity), the solution tends to 1. – Andrei Jul 09 '13 at 12:46
@Andrei, MaxSteps->Infinity – Jul 09 '13 at 12:48
@andrei, the question was not about plotting the solution, you did not provide any code for that, so I leave it up to you form the solution you need. – Jul 09 '13 at 12:50
i did in the original post, but somehow it got lost in the editing :) well nevermind that, the real problem indeed are the boundaries, because it seems like there is reflection there, and I need to simulate the system as if it were infinite – Andrei Jul 09 '13 at 12:58
2

"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... – Jens Jul 09 '13 at 18:52
@ruebenko, could you please clarify one thing? in your answer, are the boundary conditions assumed to be periodic or not? because in the documentation, for the "pseudospactral" difference order, it looks like that! I am asking because I wanted mine to be fixed.. – Andrei Jul 11 '13 at 19:14
@Andrei, Plot[vortex[xl, y] - vortex[xr, y], {y, yl, yr}] and Manipulate[ Plot[Chop[sol[xl, y, t] - vortex[xl, y] + 1], {y, yl, yr}], {t, 0, tmax}] ;-) – Jul 11 '13 at 19:36
@Andrei, but Infinity is not 5! – Jul 11 '13 at 19:37

Complex valued 2+1D nonlinear PDE using NDSolve

1 Answers1