fig8_msomi.m

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

%Note 1:  Let U=u+iv. The equations for (u, v) are
%   u_{xx}+gamma1*v_{xx}+gamma0*v+(u^2+v^2)u=mu*u
%   v_{xx}-gamma1*u_{xx}-gamma0*u+(u^2+v^2)v=mu*v
%       Operator L1 = [dxx+3*u^2+v^2-mu         gamma1*dxx+gamma0+2uv
%                      -gamma1*dxx-gamma0+2uv        dxx+3*v^2+u^2-mu]  
%   Hermitina of L1 is the transpose of L1. 

%Note 2:  In the computer execution, equations for u and v are recombined 
%into one for U=u+iv for numerical efficiency.    


Lx=40; N=128;                                     % grid 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; 
     
gamma0=0.3; gamma1=1; mu=1.2; c=1.4; DT=0.12;     
U=1.6*sech(x);                                    % initial conditions
     
for nn=1:max_iteration                            % MSOMI iterations begin
    Uold=U; 
    muold=mu;  
    L0U=(1-i*gamma1)*ifft(-kx.^2.*fft(U))-(mu+i*gamma0)*U+abs(U).^2.*U; 
    MinvL0U=ifft(fft(L0U)./(kx.^2+c)); 
    T1=real(MinvL0U); 
    T2=imag(MinvL0U);
    u=real(U); 
    v=imag(U); 
    R11=      ifft(-kx.^2.*fft(T1))+((3*u.^2+v.^2)-mu).*T1; 
    R12=-gamma1*ifft(-kx.^2.*fft(T2))+(-gamma0+2*u.*v).*T2;
    R21= gamma1*ifft(-kx.^2.*fft(T1))+( gamma0+2*u.*v).*T1;
    R22=      ifft(-kx.^2.*fft(T2))+((3*v.^2+u.^2)-mu).*T2; 
    L1HermitMinvL0U=R11+R12+i*(R21+R22); 
    MinvL1HermitMinvL0U=ifft(fft(L1HermitMinvL0U)./(kx.^2+c));
        
    if nn == 1
        U=U-MinvL1HermitMinvL0U*DT;  
        mu=mu+sum(u.*real(MinvL0U)+v.*imag(MinvL0U))*dx*DT;
    else
        T1=real(G); 
        T2=imag(G);
        F11=      ifft(-kx.^2.*fft(T1))+((3*u.^2+v.^2)-mu).*T1; 
        F12= gamma1*ifft(-kx.^2.*fft(T2))+( gamma0+2*u.*v).*T2;
        F21=-gamma1*ifft(-kx.^2.*fft(T1))+(-gamma0+2*u.*v).*T1;
        F22=      ifft(-kx.^2.*fft(T2))+((3*v.^2+u.^2)-mu).*T2; 
        L1G=F11+F12+i*(F21+F22); 
        L1GminusHU=L1G-H*U; 
        MinvL1GminusHU=ifft(fft(L1GminusHU)./(kx.^2+c)); 
        MG=ifft2(fft2(G).*(kx.^2+c)); 
        alpha=1/(sum(conj(MG).*G)*dx+abs(H).^2)-1/(sum(conj(L1GminusHU).*MinvL1GminusHU)*dx*DT); 
        theta=-sum(real(L1GminusHU).*real(MinvL0U)+imag(L1GminusHU).*imag(MinvL0U))*dx; 
        U=U+(-MinvL1HermitMinvL0U-alpha*theta*G)*DT;  
        mu=mu+(sum(u.*real(MinvL0U)+v.*imag(MinvL0U))*dx-alpha*theta*H)*DT; 
    end
    
    G=U-Uold; 
    H=mu-muold; 
    Uerror(nn)=sqrt(sum(abs(U-Uold).^2)*dx)+abs(mu-muold); Uerror(nn)
    if Uerror(nn) < error_tolerance 
	     break
	end
end
mu

figure(1)                                       % plotting results
plot(x,real(U),'b',x,imag(U),'r--', 'linewidth', 2); 
figure(2)
semilogy(1:length(Uerror), Uerror, 'linewidth', 2)