# Complex valued 2+1D nonlinear PDE using NDSolve

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.

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:

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.

Here is an example:

showStatus[status_] :=
SetNotebookStatusLine[FrontEndEvaluationNotebook[],
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]


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

• @Kuba, thanks for the edit! – user21 Jul 9 '13 at 10:05
• Your wellcome :), since You are posting a lot of differential-equations etc., You may want to look at this halirutan's script – Kuba Jul 9 '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 9 '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. – user21 Jul 9 '13 at 12:25
• "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 9 '13 at 18:52