好多研1的同学,开始用matlab进行模拟的时候都没有学过matlab,下面这个程序为光脉冲在光纤中传输时程序,其中用到了傅里叶分步法,要用傅里叶分步法的同学也可以参考借鉴一下,还有这个程序也可以用来作为锁模的机理的前序,大家参考上面的程序很容易编出锁模的程序,希望对大家有帮助,最后用的是3D输出,大家也可以参考一下3D输出的具体写法!
%===== Solving NLSE Using Split-step Fourier Method (SSFM) =======
% The NLSE is like thie:
% Uz = i/2*Uxx - i*beta*U/(1+|U|^2)
% U(z=0)= sqrt(0.11)*exp(-x.^2/0.85^2); Gauss pulse
%******************************************************************
% Split-step fast foruier method
% U(z+dz)=(L+N)U(z) L and N are operators
% L=i/2 * d/dx^2 Linear operator
% N=-i*beta/(1+|U|^2) Nonlinear operator
% Linear operator computation with fft in Fourier domain for speed
% Nonlinear operator computation in space domain
% step: 1. U-->fft(U) = ans
% 2. L(dz/2).*ans
% 3. ifft(ans)
% 4. N(dz).*ans
% 5. fft(ans)
% 6. L(dz/2).*ans
% 7. ifft(ans) =U(z+dz)
% Calculate two times L(dz/2) can get a accuracy of O(dz^3)
%******************************************************************
clear;
clc;
beta=1; gamma=0;
%******************************************************************
% Set Space & Frequency Windows
%******************************************************************
nz = 3000;
dz = 0.001;
Nstep = nz/10; % For drawing waterfall lines
xmax = 10; % x domain [-5,5]
nx = 2^9;
dx = 2*xmax/nx;
x = (-xmax:dxxmax-dx))';
fk = 2*pi/dx;
xl = length(x)/2;
k = (-xl*fk:fkxl-1)*fk)'/nx; % frequency domain [-fk/2,fk/2]
%******************************************************************
% Initialization
%******************************************************************
u0 =sech(x);
u = u0;
%******************************************************************
% Preplot
%******************************************************************
iplot = 1;
zplot(iplot) = 0;
Uplot(:,iplot) = u.*conj(u);
%******************************************************************
% Linear & Nonlinear Operators
%******************************************************************
% in fourier domain
%******************************************************************
% Split-Step Fourier Method For NLS
%******************************************************************
L=exp(i*dz/4*k.^2);
for m = 1:1:nz
%%%%%%%%%%% test for self-delfection%%%%%%%
absu2=u.*conj(u);
%u2xtemp=diff(absu2)./diff(x);
%u2x(1)=u2xtemp(1);u2x(2:nx,1)=u2xtemp(1:nx-1);
Na=i*dz*beta*u.*conj(u); % in space domain
%Nb=-Na/beta*gamma.*u2x; % in space domain
%N=exp(Na+Nb);
N=exp(Na);
%%%%%%%%%%%%%%% end test%%%%%%%
f = fftshift(fft(u));
f = L.*f;
uhalf = ifft(fftshift(f));
ua = N.*uhalf;
f = fftshift(fft(ua));
f = L.*f;
uend = ifft(fftshift(f));
% 梯形规则积分 两次迭代
u2 = exp(i/2*dz*beta*(u.*conj(u)+uend.*conj(uend))).*uhalf;
f = fftshift(fft(u2));
f = L.*f;
uend1 = ifft(fftshift(f));
u3 = exp(i/2*dz*beta.*(u.*conj(u)+uend1.*conj(uend1))).*uhalf;
f = fftshift(fft(u3));
f = L.*f;
u = ifft(fftshift(f));
er=1e-4;
u2a=uend1.*conj(uend1);
u2b=u.*conj(u);
if(abs(max(u2b-u2a))/max(u2b)) > er;
disp('Reduce step length please');break;
end
% For drawing 3-D graphic
if (rem(m,Nstep) == 0)
iplot=iplot+1;
Uplot(:,iplot)=u.*conj(u);
zplot(iplot) = m*dz;
end
end
waterfall(x,zplot,Uplot');
hidden off;
view(50,30);
title('Soliton Evolution');
xlabel('s','FontSize',18);
ylabel('\xi','FontSize',18);
zlabel('|U|^2','FontSize',14,'Rotation',0,'Position',[-xmax/1.5,-.1,max(u0)^2]);