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

计算脉冲在非线性耦合器中演化的Matlab 程序 =3w;<1 ?'  
cPNc$^Y  
%  This Matlab script file solves the coupled nonlinear Schrodinger equations of Cd"{7<OyM4  
%  soliton in 2 cores coupler. The output pulse evolution plot is shown in Fig.1 of ]2qKc  
%  Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear BR@m*JGajz  
%   pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 ceJi|`F  
usD@4!PoA  
%fid=fopen('e21.dat','w'); -dB
WpT  
N = 128;                       % Number of Fourier modes (Time domain sampling points) SnQT1U%  
M1 =3000;              % Total number of space steps 3cL iZ%6^  
J =100;                % Steps between output of space `w\P- q  
T =10;                  % length of time windows:T*T0 CCe>*tdf  
T0=0.1;                 % input pulse width fM4B.45j  
MN1=0;                 % initial value for the space output location @|c])  
dt = T/N;                      % time step )j>U4a  
n = [-N/2:1:N/2-1]';           % Index y7fy9jQ 8.  
t = n.*dt;   xi(\=LbhY  
u10=1.*sech(1*t);              % input to waveguide1 amplitude: power=u10*u10 #xw*;hW<  
u20=u10.*0.0;                  % input to waveguide 2 iI";m0Ny  
u1=u10; u2=u20;                 ~<n.5q%Z  
U1 = u1;   @|DQZt  
U2 = u2;                       % Compute initial condition; save it in U e@ZM&iR  
ww = 4*n.*n*pi*pi/T/T;         % Square of frequency. Note i^2=-1. mA+:)?e5~  
w=2*pi*n./T; ud$-A  
g=-i*ww./2;                    % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./T 3>@VPMi  
L=4;                           % length of evoluation to compare with S. Trillo's paper O.(2  
dz=L/M1;                       % space step, make sure nonlinear<0.05 SI+Uq(k  
for m1 = 1:1:M1                                    % Start space evolution ")STB8kQ  
   u1 = exp(dz*i*(abs(u1).*abs(u1))).*u1;          % 1st sSolve nonlinear part of NLS jTcv&`fAz  
   u2 = exp(dz*i*(abs(u2).*abs(u2))).*u2; -m%`Di!E  
   ca1 = fftshift(fft(u1));                        % Take Fourier transform OpEH4X.Z  
   ca2 = fftshift(fft(u2)); ()?83Xj[c  
   c2=exp(g.*dz).*(ca2+i*1*ca1.*dz);               % approximation '1gfXC  
   c1=exp(g.*dz).*(ca1+i*1*ca2.*dz);               % frequency domain phase shift   >9dD7FH  
   u2 = ifft(fftshift(c2));                        % Return to physical space lt&(S)  
   u1 = ifft(fftshift(c1)); P$#:$U@  
if rem(m1,J) == 0                                 % Save output every J steps. 1d<Uwb>  
    U1 = [U1 u1];                                  % put solutions in U array 4>>=TJ!M  
    U2=[U2 u2]; d/&> `[i  
    MN1=[MN1 m1]; '6 F-%
 
    z1=dz*MN1';                                    % output location gYc]z5`  
  end b
He' U>  
end 4}uOut  
hg=abs(U1').*abs(U1');                             % for data write to excel {4G/HW28  
ha=[z1 hg];                                        % for data write to excel 5?^
L))  
t1=[0 t']; _V-KyK  
hh=[t1' ha'];                                      % for data write to excel file 1^}I?PbqV  
%dlmwrite('aa',hh,'\t');                           % save data in the excel format LnI  
figure(1) $ItjVc@U  
waterfall(t',z1',abs(U1').*abs(U1'))               % t' is 1xn, z' is 1xm, and U1' is mxn wwB3m&  
figure(2) dWvVK("Wj  
waterfall(t',z1',abs(U2').*abs(U2'))               % t' is 1xn, z' is 1xm, and U1' is mxn gVOAB-nw  
Nhjq.&  
非线性超快脉冲耦合的数值方法的Matlab程序 W8^m-B&  
"^n,(l*4x  
在研究脉冲在非线性耦合器中的演变时，我们需要求解非线性偏微分方程组。在如下的论文中，我们提出了一种简洁的数值方法。 这里我们提供给大家用Matlab编写的计算程序。   E=p+z"Ui  
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 \:WWrY8&  
Dp ](?Yr  
o3fc-  
DVL-qt\;n  
%  This Matlab script file solves the nonlinear Schrodinger equations (?na|yd  
%  for 3 cores nonlinear coupler. The output plot is shown in Fig.2 of lb-1z]YwQ  
%  Youfa Wang and Wenfeng Wang, “A simple and effective numerical method for nonlinear 5*pzL0,Y  
%  pulse propagation in N-core optical couplers”, IEEE Photonics Technology lett. Vol.16, No.4, pp1077-1079, 2004 3S:Lce'f  
m0"K^p  
C=1;                           9Rnypzds  
M1=120,                       % integer for amplitude K -U}sW  
M3=5000;                      % integer for length of coupler -)`_w^Ox  
N = 512;                      % Number of Fourier modes (Time domain sampling points) YNEwX$)M,B  
dz =3.14159/(sqrt(2.)*C)/M3;  % length of coupler is divided into M3 segments,  make sure nonlinearity<0.05. J~k9jeq9  
T =40;                        % length of time:T*T0. l
<`>  
dt = T/N;                     % time step X g6ezl
W  
n = [-N/2:1:N/2-1]';          % Index y>P+"Z.K%}  
t = n.*dt;   I+8n;I)]X  
ww = 4*n.*n*pi*pi/T/T;        % Square of frequency. Note i^2=-1. 50^ux:Uv+N  
w=2*pi*n./T; *j%x  
g1=-i*ww./2; >X*tMhcb  
g2=-i*ww./2;                  % w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./TP=0; >.iF,[.[F<  
g3=-i*ww./2; M6)  G_-  
P1=0; j~Aq-8R=  
P2=0; h+FM?ct6}  
P3=1; f2i:I1 p("  
P=0; sS>b}u+v#!  
for m1=1:M1                 A9$x8x*Lt  
p=0.032*m1;                %input amplitude k {
*QU(  
s10=p.*sech(p.*t);         %input soliton pulse in waveguide 1 zfeT>S+  
s1=s10; {{?g%mQ6  
s20=0.*s10;                %input in waveguide 2 ci~#G[_$S  
s30=0.*s10;                %input in waveguide 3 o|kykxcq  
s2=s20; ,@`?I6nKy  
s3=s30; #).$o~1ht!  
p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s10(N,1)*s10(N,1))+abs(s10(1,1)*s10(1,1))));   {>f"&I<xw  
%energy in waveguide 1 cMw<
3u\  
p20=dt*(sum(abs(s20').*abs(s20'))-0.5*(abs(s20(N,1)*s20(N,1))+abs(s20(1,1)*s20(1,1))));   H3A$YkK [  
%energy in waveguide 2 h: ' |)O  
p30=dt*(sum(abs(s30').*abs(s30'))-0.5*(abs(s30(N,1)*s30(N,1))+abs(s30(1,1)*s30(1,1))));   f!9i6  
%energy in waveguide 3 m@td[^O-  
for m3 = 1:1:M3                                    % Start space evolution e8F]m`{_"  
   s1 = exp(dz*i*(abs(s1).*abs(s1))).*s1;          % 1st step, Solve nonlinear part of NLS ;w7mr1  
   s2 = exp(dz*i*(abs(s2).*abs(s2))).*s2; ] G&*HMtp  
   s3 = exp(dz*i*(abs(s3).*abs(s3))).*s3; 8>&@"j  
   sca1 = fftshift(fft(s1));                       % Take Fourier transform 95l)s],  
   sca2 = fftshift(fft(s2)); !_S#8"  
   sca3 = fftshift(fft(s3)); -50DGA,K6  
   sc1=exp(g1.*dz).*(sca1+i*C*sca2.*dz);           % 2nd step, frequency domain phase shift   mptFd  
   sc2=exp(g2.*dz).*(sca2+i*C*(sca1+sca3).*dz); Re$h6sh  
   sc3=exp(g3.*dz).*(sca3+i*C*sca2.*dz); bdg6B7%Q  
   s3 = ifft(fftshift(sc3)); PsC")JS  
   s2 = ifft(fftshift(sc2));                       % Return to physical space L:$4o  
   s1 = ifft(fftshift(sc1)); uUH4vUa  
end U'fP  
   p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s1(1,1)*s1(1,1)))); )9]a  
   p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1)))); M.W X&;>  
   p3=dt*(sum(abs(s3').*abs(s3'))-0.5*(abs(s3(N,1)*s3(N,1))+abs(s3(1,1)*s3(1,1)))); X3.zNHN5  
   P1=[P1 p1/p10]; wPlM= .Hq?  
   P2=[P2 p2/p10]; Hn|W3U  
   P3=[P3 p3/p10]; cHjQwl  
   P=[P p*p]; Pe`(9&iT.  
end ON :t"z5  
figure(1) aZFpt/.d  
plot(P,P1, P,P2, P,P3); o?`FjZ6;x  
6_CP?X+T  
转自：http://blog.163.com/opto_wang/

谢谢哈~！~