I am solving the nonlinear Schrodinger equation (NLSE), $$A_t+iA_{xx}+i|A|^2A=0$$ where $A$ is a complex valued function, which can be written as $A=ae^{i\theta}$ for $a,\theta$ real. Now, for particular sets of initial conditions, I find that phase dislocations occur, and hence there are singularities in the frequency, $\omega=-\theta_t$.
My question is, do phase singularities effect the accuracy of the model? Furthermore, physically, why can rapid variations in the phase occur in an equation that is derived under the WKB assumption? Does this mean the equation is not uniformly valid?
An example will be provided below.
I am using a split-step pseudo spectral method. My code is written in Fortran but I have found I can reproduce the results with a Matlab code I found online, which is provided below.
The initial condition takes the form $$A=sech(x)e^{iCx^2}$$
for $C$ a constant representing a chirp to the signal, which can lead to packet focusing due to dispersion. The evolution of the packet is shown in figure 1. 
Now, the phenomena in question occurs around $t=0.20$ when the amplitude $a$ goes to 0, and hence the phase is not well defined. This leads to singularities in $\omega$. A better picture of this is shown in figure 2.
We see that as the amplitude goes towards 0 at the front and rear of the group, i.e. the packet splits, the phase starts to become singular. Note, I've decimated $\omega$ by a factor 10 to make my point.
Note, the amplitude $a$ heuristically seems to vary slowly, but the phase clearly has rapid variations where the packet amplitude becomes small. Does this mean my results are no longer physically? Finally, and this is off topic, what sets the sign of $\omega$ at these singularities?
Note, to check the numerical accuracy of the scheme I confirm the energy is conserved to one part in $10^6$. Furthermore, there exist exact solutions to the NLSE, which I confirm stay solutions for long times with the resolution and step size used. Finally, the CFL condition for linear wave speeds ($u=1$ here) is obeyed. I have played around with the resolution and time steps, and find they give the same results.
Matlab script:
% This Matlab script file solves the nonlinear Schrodinger equation
% with the first order split-step Fourier method due to Hardin and Tappert.
% Andre Weideman, 1995
%
N=512; %number of fourier modes
dt = 0.001; %time step
M=300; % number of time steps
t=0:dt:M*dt-dt;
L=4*pi; % size of the domain
I=sqrt(-1);
% J=100;
h = L/N; % Space step
n = [-N/2:1:N/2-1]'; % Indices
x = n*h; % Grid points
%----------------------------------------
% IC: Chirped sech packet
u=sqrt(2)*sech(x).*exp(-I*x.^2);
%----------------------------------------
U = u; % Compute initial condition; save it in U
e = -4*n.*n*pi*pi/L/L; % Squares of wavenumbers.
for m = 1:1:M % Start time evolution
u = exp(dt*I*(abs(u).*abs(u))).*u; % Solve nonlinear part of NLS
c = fftshift(fft(u)); % Take Fourier transform
c = exp(dt*I*e).*c; % Advance in Fourier space
u = ifft(fftshift(c)); % Return to physical space
U = [U u];
end
