| 查看: 1201 | 回复: 0 | ||
[求助]
用分步傅里叶法求解光纤中非线性薛定谔方程
|
|
不知这里哪位虫子有用分步傅里叶法求解光纤中非线性薛定谔方程的MATLAB程序啊？我用Agrawal那本非线性光纤光学后面的程序加以改动，但是得出的图效果不理想，有没有哪位同学有程序啊？真心求交流。。。。。 附程序： %传输距离为1到4，其他参数不变，最后得出不同传输距离下频域光脉冲形状的比较 clear all; for m=1:4 %---Input Field profile distance = 10*m; %传输距离 beta2 =-1;%二阶色散参数 chirp = 0; % 初始啁啾参数 peakpower=1;%峰值功率 T0=1;%初始脉宽 r=1 ; %非线性系数 N=sqrt(r*peakpower*T0^2) ; alpha=1;%光纤损耗 %---set simulation parameters nt =1024; Tmax = 10; % FFT points and window size,设置采样点和时间窗口 step_num = round(20*distance*N^2) % No. of z steps to deltaz = distance/step_num; % step size in z dtau = (2*Tmax)/nt; % step size in tau %---tau and omega arrays tau = (-nt/2:nt/2-1)*dtau; % temporal grid omega = (pi/Tmax)*[(-nt/2:nt/2-1)]; % frequency grid %---store dispersive phase shifts to speedup code dispersion = exp((1i*0.5*beta2*omega.^2-alpha/2)*deltaz);% phase factor hhz = 1i*N^2*deltaz; % nonlinear phase factor uu =sqrt(peakpower)*exp(-0.5*tau.^2.*(1+1i*chirp)/T0^2); % gaussian temp = uu.*exp((abs(uu).^2).*hhz/2); % note hhz/2 subplot(2,1,1) if m==1 plot(fftshift(omega)/(2*pi),-log10(abs(temp).^2),'r') hold on else if m==2 plot(fftshift(omega)/(2*pi),-log10(abs(temp).^2),'y') hold on else if m==3 plot(fftshift(omega)/(2*pi),-log10(abs(temp).^2),'b') hold on else plot(fftshift(omega)/(2*pi),-log10(abs(temp).^2),'k') end end end %*********[ Beginning of MAIN Loop]*********** % scheme: 1/2N -> D -> 1/2N; first half step nonlinear for n=1:step_num f_temp = ifft(temp).*dispersion; uu = fft(f_temp); temp = uu.*exp((abs(uu).^2).*hhz/2); end uu = temp.*exp(-(abs(uu).^2).*hhz/2);% Final field temp = fftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi);%Final spectrum subplot(2,1,2) if m==1 plot(fftshift(omega)/(2*pi),log10(abs(temp)).^2,'r') hold on else if m==2 plot(fftshift(omega)/(2*pi),log10(abs(temp)).^2,'y') hold on else if m==3 plot(fftshift(omega)/(2*pi),log10(abs(temp)).^2,'b') hold on else plot(fftshift(omega)/(2*pi),log10(abs(temp)).^2,'k') end end end %***************[ End of MAIN Loop ]************** %----Plot output pulse shape and spectrum %TEMP(m,  =temp;end %for m=1:4 plot(fftshift(omega)/(2*pi),-20*(log10(abs(TEMP(m, ).^2)));% plot(fftshift(omega)/(2*pi),-20*(log10(abs(TEMP(m, ).^2)));%end |
回复此楼
点击这里搜索更多相关资源
