fig3_msom.m

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

%Note: In this case, the solution U for Eq. (134) is real, thus
%   L1= dxx+dyy+V-3U^2+mu,  where V=-V0(sin(x)^2+sin(y)^2). 
%The Hermitian of L1 is itself. 


Lx=10*pi; Ly=10*pi; N=128;                        % mesh parameters
max_iteration=1e4; 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=5; c=2.9; DT=1.7;  
U=1.15*sech(X.^2+Y.^2).*cos(X).*cos(Y);           % initial conditions
 
for nn=1:max_iteration                            % MSOM iterations start
    Uold=U; 
    L0U=ifft2(-(KX.^2+KY.^2).*fft2(U))+(V+mu).*U-U.^3; 
    MinvL0U=ifft2(fft2(L0U)./(KX.^2+KY.^2+c)); 
    L1HermitMinvL0U=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0U))+(V-3*U.^2+mu).*MinvL0U; 
    MinvL1HermitMinvL0U=ifft2(fft2(L1HermitMinvL0U)./(KX.^2+KY.^2+c));
            
    if nn == 1
        U=U-MinvL1HermitMinvL0U*DT; 
    else
        L1G=ifft2(-(KX.^2+KY.^2).*fft2(G))+(V-3*U.^2+mu).*G; 
        MinvL1G=ifft2(fft2(L1G)./(KX.^2+KY.^2+c)); 
        MG=ifft2(fft2(G).*(KX.^2+KY.^2+c)); 
        alpha=1/sum(sum(MG.*G))-1/(DT*sum(sum(L1G.*MinvL1G))); 
        U=U-(MinvL1HermitMinvL0U-alpha*sum(sum(G.*L1HermitMinvL0U))*G)*DT; 
    end
    
    G=U-Uold; 
    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
mesh(x, y, U)
figure(2)            
semilogy(1:length(Uerror), Uerror, 'linewidth', 2)