[分享]求解光孤子或超短脉冲耦合方程的Matlab程序 [复制链接]

计算脉冲在非线性耦合器中演化的Matlab 程序 `# ^0cW  
S
SmHEy*r)  
%  This Matlab script file solves the coupled nonlinear Schrodinger equations of )^/0cQcJ  
%  soliton in 2 cores coupler. The output pulse evolution plot is shown in Fig.1 of ]J@/p:S
>  
%  Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear ngUHkpYS5  
%   pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 J(iV0LAZb  
k4y}&?$B  
%fid=fopen('e21.dat','w'); `|Fp^gM  
N = 128;                       % Number of Fourier modes (Time domain sampling points) jv&+<j`r  
M1 =3000;              % Total number of space steps lhPGE_\  
J =100;                % Steps between output of space 5 9-!6;T  
T =10;                  % length of time windows:T*T0 '^}+Fv<O  
T0=0.1;                 % input pulse width (3%t+aqq  
MN1=0;                 % initial value for the space output location -cfx2;68  
dt = T/N;                      % time step +nU.p/cK+\  
n = [-N/2:1:N/2-1]';           % Index ]P1YHw9  
t = n.*dt;   oNYZIk:  
u10=1.*sech(1*t);              % input to waveguide1 amplitude: power=u10*u10 h?j_Ry  
u20=u10.*0.0;                  % input to waveguide 2 y0IK,W'&?  
u1=u10; u2=u20;                 fN[8N$1-  
U1 = u1;   !7 _\P7M  
U2 = u2;                       % Compute initial condition; save it in U #n7Yr,|Z  
ww = 4*n.*n*pi*pi/T/T;         % Square of frequency. Note i^2=-1. TV:<TR  
w=2*pi*n./T; &drFQ|  
g=-i*ww./2;                    % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./T I>n g`  
L=4;                           % length of evoluation to compare with S. Trillo's paper wE Qi0!  
dz=L/M1;                       % space step, make sure nonlinear<0.05 
V4K'R2t  
for m1 = 1:1:M1                                    % Start space evolution $>w/
Cy  
   u1 = exp(dz*i*(abs(u1).*abs(u1))).*u1;          % 1st sSolve nonlinear part of NLS Y&f\VNlT  
   u2 = exp(dz*i*(abs(u2).*abs(u2))).*u2; (tCib
4  
   ca1 = fftshift(fft(u1));                        % Take Fourier transform f/ahwz  
   ca2 = fftshift(fft(u2)); ijW7c+yd  
   c2=exp(g.*dz).*(ca2+i*1*ca1.*dz);               % approximation Lj 8<'"U#  
   c1=exp(g.*dz).*(ca1+i*1*ca2.*dz);               % frequency domain phase shift   x';uCKWV  
   u2 = ifft(fftshift(c2));                        % Return to physical space o`?zF+M0  
   u1 = ifft(fftshift(c1)); EzT`,#b  
if rem(m1,J) == 0                                 % Save output every J steps. ;l!`C':'  
    U1 = [U1 u1];                                  % put solutions in U array vsMmCd)7U  
    U2=[U2 u2]; n=!uNu7  
    MN1=[MN1 m1]; TFH&(_b  
    z1=dz*MN1';                                    % output location S`=WF^  
  end f j<H6|3  
end Ge \["`;i  
hg=abs(U1').*abs(U1');                             % for data write to excel $3;Upgv  
ha=[z1 hg];                                        % for data write to excel I/uy>*  
t1=[0 t']; !I8f#'p  
hh=[t1' ha'];                                      % for data write to excel file I1=(.
*B}  
%dlmwrite('aa',hh,'\t');                           % save data in the excel format ;YH[G;aJ  
figure(1)
qqOFr!)g  
waterfall(t',z1',abs(U1').*abs(U1'))               % t' is 1xn, z' is 1xm, and U1' is mxn 9- )qZ  
figure(2) {IM! Wb  
waterfall(t',z1',abs(U2').*abs(U2'))               % t' is 1xn, z' is 1xm, and U1' is mxn $c9k*3{<+A  
PCE4W^ns  
非线性超快脉冲耦合的数值方法的Matlab程序 J;QUPpHZ  
Pe ~c  
在研究脉冲在非线性耦合器中的演变时，我们需要求解非线性偏微分方程组。在如下的论文中，我们提出了一种简洁的数值方法。 这里我们提供给大家用Matlab编写的计算程序。   l-O$m  
Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 ls|LCQPx  
}iww:H-1  
bB6[Xj{  
Qn+:/zA;  
%  This Matlab script file solves the nonlinear Schrodinger equations EX "|H.(  
%  for 3 cores nonlinear coupler. The output plot is shown in Fig.2 of M$S]}   
%  Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear D"l+iVbBP  
%  pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 7@;">`zvm  
:1aL ?  
C=1;                           +4)7j&L  
M1=120,                       % integer for amplitude |a(fejO3  
M3=5000;                      % integer for length of coupler Fx#jV\''s  
N = 512;                      % Number of Fourier modes (Time domain sampling points) 9F##F-
%x  
dz =3.14159/(sqrt(2.)*C)/M3;  % length of coupler is divided into M3 segments,  make sure nonlinearity<0.05. -$-8W  
T =40;                        % length of time:T*T0. h*l&RR:i  
dt = T/N;                     % time step 6|;Uq'  
n = [-N/2:1:N/2-1]';          % Index \caH pof  
t = n.*dt;   GDhM<bVqM*  
ww = 4*n.*n*pi*pi/T/T;        % Square of frequency. Note i^2=-1. eSy(~Y
 
w=2*pi*n./T; )&W**!(C  
g1=-i*ww./2; L^0v\  
g2=-i*ww./2;                  % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./TP=0; p{tK_ZBy]c  
g3=-i*ww./2; p,!$/Q+l  
P1=0; >fs2kha  
P2=0; lK(Fg  
P3=1; H3KTir"on  
P=0; lj[,|[X7`  
for m1=1:M1                 c:hK$C)T  
p=0.032*m1;                %input amplitude ]k%PG-9  
s10=p.*sech(p.*t);         %input soliton pulse in waveguide 1 4<S'  
s1=s10; mY-hN|  
s20=0.*s10;                %input in waveguide 2 (?i[jO||B  
s30=0.*s10;                %input in waveguide 3 Akk 3 Qx  
s2=s20; "8<K'zeS8  
s3=s30; M"Y0jQ(  
p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s10(N,1)*s10(N,1))+abs(s10(1,1)*s10(1,1))));   0Y+FRB]u  
%energy in waveguide 1 K`6z&*  
p20=dt*(sum(abs(s20').*abs(s20'))-0.5*(abs(s20(N,1)*s20(N,1))+abs(s20(1,1)*s20(1,1))));   F:g=i}7  
%energy in waveguide 2 2xxB\J  
p30=dt*(sum(abs(s30').*abs(s30'))-0.5*(abs(s30(N,1)*s30(N,1))+abs(s30(1,1)*s30(1,1))));   0!GAk  
%energy in waveguide 3 nb,2,H  
for m3 = 1:1:M3                                    % Start space evolution F jrINxL7^  
   s1 = exp(dz*i*(abs(s1).*abs(s1))).*s1;          % 1st step, Solve nonlinear part of NLS &"E lm  
   s2 = exp(dz*i*(abs(s2).*abs(s2))).*s2; oh-|'5+,;h  
   s3 = exp(dz*i*(abs(s3).*abs(s3))).*s3; w=_Jc8/.  
   sca1 = fftshift(fft(s1));                       % Take Fourier transform i'HQQWd  
   sca2 = fftshift(fft(s2)); pV\YG B+  
   sca3 = fftshift(fft(s3)); Va<eusl  
   sc1=exp(g1.*dz).*(sca1+i*C*sca2.*dz);           % 2nd step, frequency domain phase shift   _M5%V>HO  
   sc2=exp(g2.*dz).*(sca2+i*C*(sca1+sca3).*dz); >,5i60Q  
   sc3=exp(g3.*dz).*(sca3+i*C*sca2.*dz); .qD@ Y3-  
   s3 = ifft(fftshift(sc3)); S-Fo  
   s2 = ifft(fftshift(sc2));                       % Return to physical space N/F
$bv  
   s1 = ifft(fftshift(sc1)); %V_-%/3Z  
end FY'dJY3O  
   p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s1(1,1)*s1(1,1)))); <z)m%*lvU  
   p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1)))); D]03eu  
   p3=dt*(sum(abs(s3').*abs(s3'))-0.5*(abs(s3(N,1)*s3(N,1))+abs(s3(1,1)*s3(1,1)))); .2:\:H~3  
   P1=[P1 p1/p10]; )P Jw+5  
   P2=[P2 p2/p10]; U.oksD9v  
   P3=[P3 p3/p10]; *VeW?mY,P  
   P=[P p*p]; JMa3btLy(  
end E1V^}dn  
figure(1) Mt>oI SN&d  
plot(P,P1, P,P2, P,P3); Zj9c9  
uGH?N  
转自：http://blog.163.com/opto_wang/

谢谢哈~！~