在这篇文章中，我们给出了非线性脉冲和孤子耦合的理论，它包含了2阶3阶色散，模间色散，损耗增益，色散不匹配，高阶非线性等等。文章主要说明了几个非线性耦合的基本问题： 非线性脉冲耦合行为主要依赖于参量Lcd （色散长度乘耦合系数），而不是输入脉冲形状，因而，孤子，高斯脉冲以及类高斯脉冲表现出同样的开关特性。根据Lcd, 耦合器可分为三个工作区，每个区有自己独特的耦合行为。一个脉冲是否表现为连续波的耦合行为还是超快脉冲的耦合行为，决定于Lcd 而不是脉宽。 我们显示0.8皮秒的脉冲也遵从连续波的方程和耦合行为， 可以当连续波一样来分析。

文章参见：http://www.qiji.cn/eprint/abs/3115.html 

或 Journal of Lightwave Technology, Vol. 24, Issue 6, pp. 2458- (2006)

 

%  This Matlab script file solves the nonlinear Schrodinger equation
%  study input peak power vursus output total power for different Lcd, P~P CURVES 
%  Simulate Fig.2,Fig.3 and Fig.6 Using wangs' equations (high order dispersion and IMD);
lcd=0.2;                        % 2nd dispersion length x coupling coefficient
ldd=-1300000000;                % 3rd dispersion length x coupling coefficient
lmd=-0.;                        % intermode dispersion (IMD)
M1=43,                          % number defines amplitude of input pulse
N =512;                         % Time domain sampling points
dz =0.5*0.0031415926/4;         % space step, make sure nonlinear<0.05
M3=4000;                        % dz*M3 is normalized length, M3*dz=3.14159/2 corresponds a half beat length coupler
                                % one beat length coupler: dz*M3=3.14159. usual coupler length is a half beat.
T =40;                          % length of time:T*T0, it can.t be too big ot too small, it affect accuracy
T0=60;                          % input pulse width
delt=-0.0;                      % normalized dispersion mismatch=(b1-b2)/2k
Ld=T0^2/20*1000;                % dispersion length=T0^2/abs(b2),b2=-20ps^2/Km for 1.55um
dt = T/N;                       % time step
n = [-N/2:1:N/2-1]';            % Indices
t = n.*dt; 
ww = 4*n.*n*pi*pi/T/T;          % Square of frequency. Note i^2=-1.
w=2*pi*n./T;
www=w.*ww;
g1=i*(delt/T0)*w-i*ww./(2*lcd)-i*www./(6*ldd);
g2=i*(delt/T0)*w-i*ww./(2*lcd)-i*www./(6*ldd);         %w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./TP=0;
P1=0;
P2=1;
P=0;
for m1=1:1:M1
p=0.5+0.03*m1,
%s10=p.*sech(p.*t).*exp(0.5*sqrt(-1.)*t.*t);
%s10=p.*sech(t);                                        % input
%s10=p.*exp(-0.5*(1+sqrt(-1.))*t.*t);
s10=p.*exp(-t.*t./(2*1.146*1.146));                  
%s10=p.*exp(-t.*t./2);
s1=s10;
s20=0*s10;
s2=s20;
p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s10(N,1)*s10(N,1))+abs(s10(1,1)*s10(1,1)))); %energy in waveguide 1
p20=dt*(sum(abs(s20').*abs(s20'))-0.5*(abs(s20(N,1)*s20(N,1))+abs(s20(1,1)*s20(1,1)))); %energy in waveguide 2
for m3 = 1:1:M3                                    % Start space evolution
   s1 = exp(dz*i*4*(abs(s1).*abs(s1))).*s1;        % Solve nonlinear part of NLS
   s2 = exp(dz*i*4*(abs(s2).*abs(s2))).*s2;
   sca1 = fftshift(fft(s1));                       % Take Fourier transform
   sca2 = fftshift(fft(s2));
   sc2=exp(g2.*dz).*(sca2+i*(1+lmd*w).*sca1.*dz);
   sc1=exp(g1.*dz).*(sca1+i*(1+lmd*w).*sca2.*dz);  % frequency domain phase shift   % Advance in Fourier space
   s2 = ifft(fftshift(sc2));                       % Return to physical space
   s1 = ifft(fftshift(sc1));
 end
   p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s1(1,1)*s1(1,1))));
   p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1))));
   P1=[P1 p1/p10];
   P2=[P2 p2/p10];
   P=[P p*p];
end
figure(1)
plot(P,P1, P,P2);