好多研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]);