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

计算脉冲在非线性耦合器中演化的Matlab 程序 )NW6?Pu"  
2j&@p>  
%  This Matlab script file solves the coupled nonlinear Schrodinger equations of ~oD8Rnf  
%  soliton in 2 cores coupler. The output pulse evolution plot is shown in Fig.1 of
)@g;j>  
%  Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear `$5UH
a2/  
%   pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 zA+@FR?  
-=}3j&,\R  
%fid=fopen('e21.dat','w'); tpf7_YP_!-  
N = 128;                       % Number of Fourier modes (Time domain sampling points) g:)DNy  
M1 =3000;              % Total number of space steps 1(dj[3Mt  
J =100;                % Steps between output of space Oe]&(  
T =10;                  % length of time windows:T*T0 iU2KEqCm  
T0=0.1;                 % input pulse width ~=n#}{/  
MN1=0;                 % initial value for the space output location %!j:fJ()  
dt = T/N;                      % time step @ GDX7TPV  
n = [-N/2:1:N/2-1]';           % Index @e2}BhB2  
t = n.*dt;   viaJblYj(f  
u10=1.*sech(1*t);              % input to waveguide1 amplitude: power=u10*u10 9}tG\0tL*  
u20=u10.*0.0;                  % input to waveguide 2 \ZXLX'-  
u1=u10; u2=u20;                 'ktHPn ,K  
U1 = u1;   2 YxTMT  
U2 = u2;                       % Compute initial condition; save it in U `k{& /]  
ww = 4*n.*n*pi*pi/T/T;         % Square of frequency. Note i^2=-1. x;E2~&E  
w=2*pi*n./T; :
os
z  
g=-i*ww./2;                    % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./T QBJ3iQs1  
L=4;                           % length of evoluation to compare with S. Trillo's paper quu*xJ;Ci  
dz=L/M1;                       % space step, make sure nonlinear<0.05 c'fSu;1  
for m1 = 1:1:M1                                    % Start space evolution 90N`CXas
 
   u1 = exp(dz*i*(abs(u1).*abs(u1))).*u1;          % 1st sSolve nonlinear part of NLS %"$@%"8;3  
   u2 = exp(dz*i*(abs(u2).*abs(u2))).*u2; l5t2\Fl  
   ca1 = fftshift(fft(u1));                        % Take Fourier transform $ChK]v 6C  
   ca2 = fftshift(fft(u2)); C*6S@4k  
   c2=exp(g.*dz).*(ca2+i*1*ca1.*dz);               % approximation u' Qd,  
   c1=exp(g.*dz).*(ca1+i*1*ca2.*dz);               % frequency domain phase shift   vwGeD|Fb5  
   u2 = ifft(fftshift(c2));                        % Return to physical space E8}+k o  
   u1 = ifft(fftshift(c1)); C'mL&  
if rem(m1,J) == 0                                 % Save output every J steps. #VbVsl  
    U1 = [U1 u1];                                  % put solutions in U array 0Fr1Ku!  
    U2=[U2 u2]; ,d,\-x-+/  
    MN1=[MN1 m1]; !>^JSHR4t  
    z1=dz*MN1';                                    % output location Wa"(m*hW  
  end HL{$ ^l#v  
end hq>Csj==@  
hg=abs(U1').*abs(U1');                             % for data write to excel 
V_^@  
ha=[z1 hg];                                        % for data write to excel Z'v-F^  
t1=[0 t']; mryN}  
hh=[t1' ha'];                                      % for data write to excel file kAzd8nJ'  
%dlmwrite('aa',hh,'\t');                           % save data in the excel format d.<~&.-$  
figure(1) 4/>Our 5  
waterfall(t',z1',abs(U1').*abs(U1'))               % t' is 1xn, z' is 1xm, and U1' is mxn 748CD{KxW  
figure(2) h rN%  
waterfall(t',z1',abs(U2').*abs(U2'))               % t' is 1xn, z' is 1xm, and U1' is mxn w=b(X q+:  
2h^WYpCm  
非线性超快脉冲耦合的数值方法的Matlab程序 ,t$,idcT+
 
JN3cg  
在研究脉冲在非线性耦合器中的演变时，我们需要求解非线性偏微分方程组。在如下的论文中，我们提出了一种简洁的数值方法。 这里我们提供给大家用Matlab编写的计算程序。   uV6g[J  
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 ,2[ra9n  
Yn51U6_S  
ffDc6*.Q  
i^z`"3#LE  
%  This Matlab script file solves the nonlinear Schrodinger equations !mfJpJ  
%  for 3 cores nonlinear coupler. The output plot is shown in Fig.2 of ,\PVC@xJ  
%  Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear Zy"=y+e!E;  
%  pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 MFit
|C  
0(>rG{u  
C=1;                           6iezLG5  
M1=120,                       % integer for amplitude Bn wzcl  
M3=5000;                      % integer for length of coupler h+7>#*DH  
N = 512;                      % Number of Fourier modes (Time domain sampling points) eOlKbJU  
dz =3.14159/(sqrt(2.)*C)/M3;  % length of coupler is divided into M3 segments,  make sure nonlinearity<0.05. .5HQ   
T =40;                        % length of time:T*T0. Al"3 kRJJ  
dt = T/N;                     % time step :Waox"#=g  
n = [-N/2:1:N/2-1]';          % Index 9|r* pK[  
t = n.*dt;   Ps[$.h  
ww = 4*n.*n*pi*pi/T/T;        % Square of frequency. Note i^2=-1. +RZ~LA\+  
w=2*pi*n./T; yf1CXldi  
g1=-i*ww./2; V-O(U*]  
g2=-i*ww./2;                  % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./TP=0; VkmRh,T  
g3=-i*ww./2; g;p)n  
P1=0; MEZ{j%-a  
P2=0; KlxN~/gyik  
P3=1; | FM }  
P=0; #} ,x @]
p  
for m1=1:M1                 3-Bl  
p=0.032*m1;                %input amplitude aC=['a>)  
s10=p.*sech(p.*t);         %input soliton pulse in waveguide 1 |(IO=V4P  
s1=s10; Q%ad q-B  
s20=0.*s10;                %input in waveguide 2 'JmBh@A  
s30=0.*s10;                %input in waveguide 3 ?2J?XS>  
s2=s20; T`YwJ6N  
s3=s30; Jn}n*t3  
p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s10(N,1)*s10(N,1))+abs(s10(1,1)*s10(1,1))));   0NE{8O0;Fr  
%energy in waveguide 1 pgc3jP!  
p20=dt*(sum(abs(s20').*abs(s20'))-0.5*(abs(s20(N,1)*s20(N,1))+abs(s20(1,1)*s20(1,1))));   vn').\,P2O  
%energy in waveguide 2 U..<iNQE5  
p30=dt*(sum(abs(s30').*abs(s30'))-0.5*(abs(s30(N,1)*s30(N,1))+abs(s30(1,1)*s30(1,1))));   h^|5|l  
%energy in waveguide 3 'A{h iY  
for m3 = 1:1:M3                                    % Start space evolution =jAFgwP\  
   s1 = exp(dz*i*(abs(s1).*abs(s1))).*s1;          % 1st step, Solve nonlinear part of NLS w_-+o^  
   s2 = exp(dz*i*(abs(s2).*abs(s2))).*s2; X~U >LLr  
   s3 = exp(dz*i*(abs(s3).*abs(s3))).*s3; m
O rWJ~=  
   sca1 = fftshift(fft(s1));                       % Take Fourier transform #B}?Zg  
   sca2 = fftshift(fft(s2)); I+?hG6NM  
   sca3 = fftshift(fft(s3)); _]>JB0IY  
   sc1=exp(g1.*dz).*(sca1+i*C*sca2.*dz);           % 2nd step, frequency domain phase shift   C*~aSl7  
   sc2=exp(g2.*dz).*(sca2+i*C*(sca1+sca3).*dz); %IZ)3x3l  
   sc3=exp(g3.*dz).*(sca3+i*C*sca2.*dz); !>.vh]8g  
   s3 = ifft(fftshift(sc3)); M].8HwC+  
   s2 = ifft(fftshift(sc2));                       % Return to physical space 9(1rh9`=  
   s1 = ifft(fftshift(sc1)); OKue" p  
end !XE aF]8  
   p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s1(1,1)*s1(1,1)))); iw]k5<qKj  
   p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1)))); +c,[ Q  
   p3=dt*(sum(abs(s3').*abs(s3'))-0.5*(abs(s3(N,1)*s3(N,1))+abs(s3(1,1)*s3(1,1)))); HxwlYx,4  
   P1=[P1 p1/p10]; HOW7cV'X  
   P2=[P2 p2/p10]; fv'4f$U  
   P3=[P3 p3/p10]; fib#CY  
   P=[P p*p]; Utl t<  
end ?m%h`<wgMc  
figure(1) ISqfU]>[  
plot(P,P1, P,P2, P,P3); 19
u =W(  
J1F{v)T'?  
转自：http://blog.163.com/opto_wang/

谢谢哈~！~