matlab 非线性仿真,科学网—含高阶色散和高阶非线性项的非线性耦合仿真Matlab源程序 - 王又法的博文...

如下Matlab 程序计算 论文Youfa Wang and Wenfeng Wang,  Study of ultrafast pulse coupling dynamics considering retarded

nonlinear response and self-steepening effects, Joutnal of lightwave Technology, Vol.24, No.2,

pp1041-1047,中的 图2。 做一点简单改动，文中其他图，也可以计算。

% nonlinear equation including nonlinear delay, self-steepening, and higher order dispersion

%  Youfa Wang and Wenfeng Wang,  Study of ultrafast pulse coupling dynamics considering retarded

%  nonlinear response and self-steepening effects, Joutnal of lightwave Technology, Vol.24, No.2,

%  pp1041-1047,

MN3=0;

lcd=0.5;                    %lcd>0 corresponding to beita2<0, lcd<0 corresponds

%lcd=1000;

ldd=-200*lcd;

lmd=-0.00;

fr=0.18;

st=0.028;                         %st=1/(w0*t0)

%st=0;

J=400;

N =1024;                          % Number of Fourier modes (Time domain sampling points)

M3=6000;

dz =1.570796326/M3;               % space step, make sure nonlinear<0.05

T =20;                            % normalized time, t=T*T0, it can't be too big ot too small, it affect accuracy

delt=-0.0;                        % normal mismatch=(b1-b2)/2c

T0=28.4;

dt = T/N;                         % time step

n = [-N/2:1:N/2-1]';              % Indices

t = n.*dt;

%----Ranman response in fiber---------------------------

n1=[0:1:N/2]';         % the response is 0, if t1<0. it is not easy to use t. the result in funcation no difference

t1=n1*dt;

ht=(32*32+12.2*12.2)/(32*32*12.2)*T0*exp(-t1*T0/32).*sin(t1*T0/12.2);

%----------------display convolution-------------

y=[-N/2:1:N-1]';

t2=dt*y;

%--------------------------------------------------------------

ww = 4*n.*n*pi*pi/T/T;            % Square of frequency. Note i^=-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;

p=1.2;

s10=p.*sech(t);

s1=s10;

s20=0*s10;

s2=s20;

S1=s1.*0; S2=s1.*0;SC1=s1.*0; SC2=s1.*0;

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

for m3 = 1:1:M3                                    % Start space evolution

s1squre=abs(s1.*s1);

s1cubic=s1squre.*s1;

s2squre=abs(s2.*s2);

s2cubic=s2squre.*s2;

chta=conv(ht,s1squre);

cht1=dt*chta(1:N);

chtb=conv(ht,s2squre);

cht2=dt*chtb(1:N);

s1(1)=s1(1)+4*(1-fr)*i*dz*s1cubic(1)-4*(1-fr)*st*dz/dt*s1cubic(1)+4*fr*i*dz*s1(1).*cht1(1)-...

4*fr*st/dz*s1(1).*cht1(1);

s1(2:N)=s1(2:N)+4*(1-fr)*i*dz*s1cubic(2:N)-4*(1-fr)*st*dz/dt*(s1cubic(2:N)-...

s1cubic(1:N-1))+4*fr*i*dz*s1(2:N).*cht1(2:N)-4*fr*st*dz/dt*(s1(2:N).*cht1(2:N)-s1(1:N-1).*cht1(1:N-1));

s2(1)=s2(1)+4*(1-fr)*i*dz*s2cubic(1)-4*(1-fr)*st*dz/dt*s2cubic(1)+4*fr*i*dz*s2(1).*cht2(1)-...

4*fr*st/dz*s2(1).*cht2(1);

s2(2:N)=s2(2:N)+4*(1-fr)*i*dz*s2cubic(2:N)-4*(1-fr)*st*dz/dt*(s2cubic(2:N)-...

s2cubic(1:N-1))+4*fr*i*dz*s2(2:N).*cht2(2:N)-4*fr*st*dz/dt*(s2(2:N).*cht2(2:N)-s2(1:N-1).*cht2(1:N-1));

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

s2 = ifft(fftshift(sc2));                       % Return to physical space

s1 = ifft(fftshift(sc1));

% if rem(m3,J) == 0                              % Save output every J steps.

% S1 = [S1 s1];                                  % put solutions in U array

% S2=[S2 s2];

% SC1 = [SC1 sc1];                               % put solutions in U array

% SC2=[SC2 sc2];

% MN3=[MN3 m3];

% z3=dz*MN3';                                    % output location

%end

m3,a=max(abs(cht1)),

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+p2-p10,

%figure(3)

%waterfall(w',z3',abs(SC1').*abs(SC1'))  %t' is 1xn, z' is 1xm, and U1' is mxn

%figure(4)

%waterfall(w',z3',abs(SC2').*abs(SC2'))  %t' is 1xn, z' is 1xm, and U1' is mxn

figure(1)

plot(t,abs(s1'.*s1'), t,abs(s2'.*s2'),'.');

figure(2)

plot(t,abs(s2'.*s2'),'.');

%figure(2)

%plot(w,abs(sc1'.*sc1'), w,abs(sc2'.*sc2'));

%-----------display convolution---

%plot(t2,chta)

%------------------------------------