fig2_som.m

用SOM方法计算非线性shrodinger方程

%Note 1:  Let U=u+iv. The equations for (u, v) are Eqs. (135)-(136) in the paper. 
%       Operator L1 = [dxx+dyy+V+3u^2+v^2+mu           2uv
%                             2uv              dxx+dyy+V+3v^2+u^2+mu],
%       where V=-V0(sin(x)^2+sin(y)^2). 
%   Hermitina of L1 is L1 itself. 

%Note 2:  In the computer execution, equations for u and v are recombined 
%into one for U=u+iv for numerical efficiency. In doing so, for F=f+ig, 
% L1[f, g]' is recombined into  Fxx+Fyy+(V+|U|^2+mu)F+2U*real(conj(U)*F)


Lx=12*pi; Ly=12*pi; N=256;                          % mesh parameters
max_iteration=1e5; error_tolerance=1e-10;        
x=-Lx/2:Lx/N:Lx/2-Lx/N; dx=Lx/N; kx=[0:N/2-1  -N/2:-1]*2*pi/Lx;
y=-Ly/2:Ly/N:Ly/2-Ly/N; dy=Ly/N; ky=[0:N/2-1  -N/2:-1]*2*pi/Ly;
[X, Y]=meshgrid(x, y); [KX, KY]=meshgrid(kx, ky);  

V=-6*(sin(X).^2+sin(Y).^2); mu=3; c=3.7; DT=0.8;    % other parameters and initial conditions   
U=1.7*(exp(-X.^2-Y.^2)+i*exp(-(X-pi).^2-Y.^2)-exp(-(X-pi).^2-(Y-pi).^2)-i*exp(-X.^2-(Y-pi).^2)); 
                      
for nn=1:max_iteration                              % SOM iterations start 
    Uold=U; 
    L0U=ifft2(-(KX.^2+KY.^2).*fft2(U))+(V+abs(U).^2+mu).*U; 
    MinvL0U=ifft2(fft2(L0U)./(KX.^2+KY.^2+c)); 
    L1HermitMinvL0U=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0U))+(V+abs(U).^2+mu).*MinvL0U+2*U.*real(conj(U).*MinvL0U); 
    MinvL1HermitMinvL0U=ifft2(fft2(L1HermitMinvL0U)./(KX.^2+KY.^2+c));
    U=U-MinvL1HermitMinvL0U*DT; 
    Uerror(nn)=sqrt(sum(sum(abs(U-Uold).^2))*dx*dy); Uerror(nn)
    if Uerror(nn) < error_tolerance 
         break
    end
end

figure(1)                                           % plotting results
imagesc(x, y, abs(U))
axis([-2.5*pi 3.5*pi -2.5*pi 3.5*pi])
figure(2)
imagesc(x, y, angle(U))
axis([-2.5*pi 3.5*pi -2.5*pi 3.5*pi])
figure(3)            
semilogy(1:length(Uerror), Uerror, 'linewidth', 2)