fig1.m

来自「计算非线性shrodinger方程」· M 代码 · 共 33 行

Lx=10*pi; Ly=10*pi; N=128;                        % mesh parameters
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); K2=KX.^2+KY.^2;  

mu=3.7; Dt=1;                              
u=1.5*exp(-(X.^2+Y.^2));                          % initial condition
uerror=1;  Uerror=[];                             % initialize E_n

while uerror >= 1e-10
    u_old=u; DEL_u=real(-ifft2(K2.*fft2(u)));
    if uerror >= 10^(-3)        % when E_n > 10^(-3), compute c and gamma
        dVu=2*u.^3;   u_u=sum(sum(u.^2));
        u_Lu=sum(sum(dVu.*u));   DELu_DELu=sum(sum(DEL_u.^2));
        DELu_Lu=sum(sum(dVu.*DEL_u));   u_DELu=sum(sum(u.*DEL_u));
        c=(u_Lu*DELu_DELu-DELu_Lu*u_DELu)/(u_Lu*u_DELu-DELu_Lu*u_u); 
        u_Nu=c*u_u-u_DELu;   alpha=u_Lu/u_Nu;   gamma=1+1/(alpha*Dt);
      else       % once E_n < 10^(-3), use previously computed c and gamma.
        u_Nu=sum(sum(c*u.^2-u.*DEL_u));
    end
    L0u=DEL_u + (3*(cos(X).^2+cos(Y).^2)-mu).*u + u.^3;
    u=u+Dt*real(ifft2(fft2(L0u)./(c+K2))-u*gamma*sum(sum(u.*L0u))/u_Nu);
    uerror=sqrt(sum(sum((u-u_old).^2))*dx*dy)     % new E_n
    Uerror=[Uerror uerror]; 
end

figure(1)                                         % plotting results
mesh(x, y, u)
figure(2)            
semilogy(1:length(Uerror), Uerror, 'linewidth', 2)