分布傅里叶方法求解光纤中的非线性薛定谔方程仿真

% This code solves the NLS equation with the split-step method %本代码使用分布傅里叶方法解决了光纤中的NLS非线性薛定谔方程 % idu/dz - sgn(beta2)/2 d^2u/d(tau)^2 + N^2*|u|^2*u=0 % 此为归一化NLS，N是孤子阶数，sgn是阶跃函数，u是慢变振幅，beta2是色散系数，tau是时间T/脉冲宽度 % Written by Govind P. Agrawal in March 2005 for the NLFO book %---Specify input parameters clear all; distance = 10;('Enter fiber length (in units of L_D) ='); %以色散长度LD为单位的光纤长度 beta2 = -1;('dispersion: 1 for normal, -1 for anomalous'); %1为正常色散，-1为反常色散 N = 2 ;('Nonlinear parameter N = '); % Soliton order 孤子阶数 mshape = 0 ;('m = 0 for sech, m > 0 for super-Gaussian = ');%0为双曲正割形，大于0为超高斯形 chirp0 = 0; % input pulse chirp (default value) %---set simulation parameters nt = 1024; Tmax = 32; % FFT points and window size step_num = round(20*distance*N^2); % No. of z steps z方向步数是20*长度*孤子阶数平方 800步 deltaz = distance/step_num; % step size in z z步长是长度/步数 0.0125 dtau = (2*Tmax)/nt; % step size in tau 时间t步长0.0625 %---tau and omega arrays tau = (-nt/2:nt/2-1)*dtau; % temporal grid -32到32内计算的时间点t1,t2.... omega = (pi/Tmax) * [(0:nt/2-1) (-nt/2:-1)]; % frequency grid 频率间隔0-50.167 -50.265--0.0982 %---Input Field profile if mshape==0 % soliton sech shape 如果是双曲正割形 uu = sech(tau).*exp(-0.5i* chirp0* tau.^2); % else % super-Gaussian uu = exp(-0.5*(1+1i*chirp0).*tau.^(2*mshape)); end %---Plot input pulse shape and spectrum temp = fftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi); % spectrum 得出输入光谱 figure; subplot(2,1,1); plot (tau, abs(uu).^2, '--r'); hold on; %abs取模后平方 axis([-5 5 0 inf]); xlabel('Normalized Time'); ylabel('Normalized Power'); title('Input and Output Pulse Shape and Spectrum'); subplot(2,1,2); plot (fftshift(omega)/(2*pi), abs(temp).^2, '--r'); hold on; axis([-.5 .5 0 inf]); xlabel('Normalized Frequency'); ylabel('Spectral Power'); %---Store dispersive phase shifts to speedup code dispersion = exp(0.5i*beta2*omega.^2*deltaz); % phase factor 色散项，简化了三阶色散和吸收损耗，实为傅里叶变换后的exp（hD） hhz = 1i*N^2*deltaz; % nonlinear phase factor 非线性效应项，实则为exp(hN)=exp(hhz*|u|^2) % ********* [ Beginning of MAIN Loop] *********** % scheme: 1/2N -> D -> 1/2N; first half step nonlinear temp = uu.*exp(abs(uu).^2.*hhz/2); % note hhz/2 先对初始函数作用非线性N/2算符 for n=1:step_num f_temp = ifft(temp).*dispersion;%取傅里叶变换后再作用色散算符 uu = fft(f_temp);%作用完色散到达h/2处，傅里叶变换回去 temp = uu.*exp(abs(uu).^2.*hhz);%再作用非线性N算符 end uu = temp.*exp(-abs(uu).^2.*hhz/2); % Final field 最后作用非线性N/2算符，此为最终时域谱 temp = fftshift(ifft(uu)).* (nt*dtau)/sqrt(2*pi); % Final spectrum 此为最终频域谱 % *************** [ End of MAIN Loop ] ************** %----Plot output pulse shape and spectrum subplot(2,1,1) plot (tau, abs(uu).^2, '-b') subplot(2,1,2) plot(fftshift(omega)/(2*pi), abs(temp).^2, '-b')