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

计算脉冲在非线性耦合器中演化的Matlab 程序 <dD)>Y.  
Uhw:XV@m  
%  This Matlab script file solves the coupled nonlinear Schrodinger equations of /<R[X>]<F  
%  soliton in 2 cores coupler. The output pulse evolution plot is shown in Fig.1 of Irc(5rD7  
%  Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear u_Wftb?9  
%   pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 *el~sor;S  
t@;r~Sb  
%fid=fopen('e21.dat','w'); yrF"`/zv6|  
N = 128;                       % Number of Fourier modes (Time domain sampling points) ;
4'pucq5/  
M1 =3000;              % Total number of space steps m]?C @ina  
J =100;                % Steps between output of space W"v"mjYud  
T =10;                  % length of time windows:T*T0 sGp]jqX2,m  
T0=0.1;                 % input pulse width lEYAq'=  
MN1=0;                 % initial value for the space output location W{is2s  
dt = T/N;                      % time step +F4SU(T  
n = [-N/2:1:N/2-1]';           % Index 6Mj
(B*c  
t = n.*dt;   0
! W$Cz[  
u10=1.*sech(1*t);              % input to waveguide1 amplitude: power=u10*u10 (M-Wea!q  
u20=u10.*0.0;                  % input to waveguide 2 tW -f_0a.  
u1=u10; u2=u20;                 ?'IY0^  
U1 = u1;   Q H57[Yg  
U2 = u2;                       % Compute initial condition; save it in U fEB&)mM  
ww = 4*n.*n*pi*pi/T/T;         % Square of frequency. Note i^2=-1. fZtuP1-4  
w=2*pi*n./T; 1EemVZdY  
g=-i*ww./2;                    % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./T 1!
=^mu8  
L=4;                           % length of evoluation to compare with S. Trillo's paper q2e=(]rKE{  
dz=L/M1;                       % space step, make sure nonlinear<0.05 K(3_1*
e  
for m1 = 1:1:M1                                    % Start space evolution *DcB?8%  
   u1 = exp(dz*i*(abs(u1).*abs(u1))).*u1;          % 1st sSolve nonlinear part of NLS di4>Ir~]  
   u2 = exp(dz*i*(abs(u2).*abs(u2))).*u2; v;o/M6GL5  
   ca1 = fftshift(fft(u1));                        % Take Fourier transform
f.G"[p  
   ca2 = fftshift(fft(u2)); =#>F' A  
   c2=exp(g.*dz).*(ca2+i*1*ca1.*dz);               % approximation )>]@@Trx  
   c1=exp(g.*dz).*(ca1+i*1*ca2.*dz);               % frequency domain phase shift   O>3f*Cc  
   u2 = ifft(fftshift(c2));                        % Return to physical space ,<t)aZL,A;  
   u1 = ifft(fftshift(c1)); [vTk*#Cl4  
if rem(m1,J) == 0                                 % Save output every J steps. I/hq8v~
S  
    U1 = [U1 u1];                                  % put solutions in U array ms{iQ:'9  
    U2=[U2 u2]; *hIjVKTu79  
    MN1=[MN1 m1]; skP'- ^F~  
    z1=dz*MN1';                                    % output location b[rVr J  
  end C0}@0c  
end H7{I[>:  
hg=abs(U1').*abs(U1');                             % for data write to excel ZZ324UuATX  
ha=[z1 hg];                                        % for data write to excel sW&5Mu-  
t1=[0 t']; B2^*Sr[  
hh=[t1' ha'];                                      % for data write to excel file #GuN.`__n,  
%dlmwrite('aa',hh,'\t');                           % save data in the excel format z(n Ba]^[F  
figure(1) uZml.#@4  
waterfall(t',z1',abs(U1').*abs(U1'))               % t' is 1xn, z' is 1xm, and U1' is mxn 80%L!x|  
figure(2) P
47x-;  
waterfall(t',z1',abs(U2').*abs(U2'))               % t' is 1xn, z' is 1xm, and U1' is mxn >/+R~ n  
isU4D  
非线性超快脉冲耦合的数值方法的Matlab程序 NHQi_U  
ez14f$cJ+  
在研究脉冲在非线性耦合器中的演变时，我们需要求解非线性偏微分方程组。在如下的论文中，我们提出了一种简洁的数值方法。 这里我们提供给大家用Matlab编写的计算程序。   )wNcz~ Y  
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 dT7!+)s5-  
[.'9Sw  
rlQ=rNrG&E  

KA? J:  
%  This Matlab script file solves the nonlinear Schrodinger equations uqvS  
%  for 3 cores nonlinear coupler. The output plot is shown in Fig.2 of *t300`x
  
%  Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear \IP 9EFA  
%  pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 \YPvpUg  
(9Of,2]&E  
C=1;                           QTospHf`  
M1=120,                       % integer for amplitude uK=)65]  
M3=5000;                      % integer for length of coupler mRIBE9K+&  
N = 512;                      % Number of Fourier modes (Time domain sampling points) r1BL?&X-  
dz =3.14159/(sqrt(2.)*C)/M3;  % length of coupler is divided into M3 segments,  make sure nonlinearity<0.05. J,*+Ak ~  
T =40;                        % length of time:T*T0. 8?LHYdJ  
dt = T/N;                     % time step n.=Zw2F
E  
n = [-N/2:1:N/2-1]';          % Index 3}lIY7O  
t = n.*dt;   8`z  
ww = 4*n.*n*pi*pi/T/T;        % Square of frequency. Note i^2=-1. 6xLQ  
w=2*pi*n./T; >BZ,g!N,J}  
g1=-i*ww./2; L:\>)6]Ls  
g2=-i*ww./2;                  % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./TP=0; <DN7  
g3=-i*ww./2; 3<>DDY2bl  
P1=0; .q<5O
E(f  
P2=0; yRR[M@Y  
P3=1; p$}/~5b}4  
P=0; t=fr`|!  
for m1=1:M1                 _%u t#  
p=0.032*m1;                %input amplitude "hnvND4=  
s10=p.*sech(p.*t);         %input soliton pulse in waveguide 1 n.XgGT=L  
s1=s10; {_4`0J`3  
s20=0.*s10;                %input in waveguide 2
^Rh}[  
s30=0.*s10;                %input in waveguide 3 G
km{b[  
s2=s20; [)?yH3  
s3=s30; %c@PTpAM  
p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s10(N,1)*s10(N,1))+abs(s10(1,1)*s10(1,1))));   Q^8/"aV\  
%energy in waveguide 1 =E62N7_`=  
p20=dt*(sum(abs(s20').*abs(s20'))-0.5*(abs(s20(N,1)*s20(N,1))+abs(s20(1,1)*s20(1,1))));   %-[*G;c'w  
%energy in waveguide 2 B'I_i$g4w  
p30=dt*(sum(abs(s30').*abs(s30'))-0.5*(abs(s30(N,1)*s30(N,1))+abs(s30(1,1)*s30(1,1))));   _ glB<r$  
%energy in waveguide 3 LkWY6 ?$U  
for m3 = 1:1:M3                                    % Start space evolution ~gaWZQXyu  
   s1 = exp(dz*i*(abs(s1).*abs(s1))).*s1;          % 1st step, Solve nonlinear part of NLS s Fx0  
   s2 = exp(dz*i*(abs(s2).*abs(s2))).*s2; |rI;OvZ\  
   s3 = exp(dz*i*(abs(s3).*abs(s3))).*s3; OAaLCpRp  
   sca1 = fftshift(fft(s1));                       % Take Fourier transform Sx1|Oq]  
   sca2 = fftshift(fft(s2)); 1DlXsup&?#  
   sca3 = fftshift(fft(s3)); <cO `jK  
   sc1=exp(g1.*dz).*(sca1+i*C*sca2.*dz);           % 2nd step, frequency domain phase shift   ??f,(om  
   sc2=exp(g2.*dz).*(sca2+i*C*(sca1+sca3).*dz); ^VEaOKMr  
   sc3=exp(g3.*dz).*(sca3+i*C*sca2.*dz); b&6lu4D  
   s3 = ifft(fftshift(sc3)); Uy|Tu~  
   s2 = ifft(fftshift(sc2));                       % Return to physical space PZVH=dagq  
   s1 = ifft(fftshift(sc1)); M+I9k;N6&  
end 5,+fM6^V  
   p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s1(1,1)*s1(1,1)))); s4j]kH  
   p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1)))); y>cLG5v  
   p3=dt*(sum(abs(s3').*abs(s3'))-0.5*(abs(s3(N,1)*s3(N,1))+abs(s3(1,1)*s3(1,1)))); 3WpQzuHPT  
   P1=[P1 p1/p10]; )q`.tsR>  
   P2=[P2 p2/p10]; tt|P-p
-  
   P3=[P3 p3/p10]; 97/4J  
   P=[P p*p]; v\E6N2.S  
end 4"UH~A;^  
figure(1) b~2LD3"3  
plot(P,P1, P,P2, P,P3); +BDW1%  
{<1uV']x  
转自：http://blog.163.com/opto_wang/

谢谢哈~！~