分步傅里叶算法的MATLAB程序实现举例
模型：
线性部分： 非线性部分：?U1?2i?z?U2?x2?nU?0 n?U?x,z?2?d?z??2n?x2其中：
d?z??exp??z?或d?z??exp???z?
?U?z=i?2U2?x2 两边同时对x变量作傅里叶变换
?U~i2~?z=2?i??U 两边积分
?x?h21~x?h2ixdU=U~?x2?i??2dz 即
~lnU?x?h/2?U~?x?=i2?i??2*h2 最后有
U~?x?h/2?=U~?x?exp??i2h??2?i??*2??
再对x变量作作傅里叶逆变换
U?x?h/2?=FFT-1???i2h???FFT??U?x???exp??2?i??*2????
?U?z?inU 两边积分
?x?h1xxUdU=??hxindz
当h?0时
ln最后有
U?x?h??inh
U?x?U?x?h??U?x?exp?inh?
折射率部分： n?U?x,z?两边同时对x变量作傅里叶变换
22n=FFT?U??d?z?*?i??n
??2FFT?U???
n=21??d?z?~2?2n?d?z?2
?x~~再对x变量作作傅里叶逆变换
?FFT?U2??????
n=FFT-1??21??dz??????MATLAB程序实现：
clear all
delta=1; x0=1;
%%------------------- n=2048; hx=0.06;
x=(-n/2:n/2-1)*hx; hw=2*pi/(n*hx);
w=fftshift((-n/2:n/2-1)*hw); %%-------------------
q=exp(-1*(x-x0).^2/2)+ exp(-1*(x+x0).^2/2); % q=sech(x);
u1(:,1)=(abs(q).^2)'; %------------------- L=500; nm=L*100; h=L/nm;
%------------------- for j=1:nm
j;
Dz=exp(delta*j*h)
D=exp(i*((i*w).^2/2)*h/2); qstep1=ifft(D.*fft(q));
n_index=ifft(fft(abs(qstep1).^2)./(Dz*w.^2+1)); N=exp(i*n_index*h); qstep2=N.*qstep1; q=ifft(D.*fft(qstep2)); u=abs(q);
r=floor(2+(j-1)/L); u1(:,r)=u'; end
z=0:L*h:L; figure(1)
mesh(x,z,u1'); view(0,90) figure(2)
plot(x,u1(:,end),'r',x,V,'b')
虚时间变换：
作虚时间变换：
z??iz
得到
?U1?2U???nU?pR(x)U?0?z2?x2 22?nn?U?x,z??d2?xMATLAB程序实现：
?U1?2U=线性部分： 2?z2?x两边同时傅里叶变换
~?U12=?i??U ?z2~两边积分
?即
x?h2xdU=?~xU1~x?h212?i??dz 2U?x?h/2?1h2ln=i?* ??~22Ux~??~最后有
h?2?1U?x?h/2?=U?x?exp??i??*?
2??2~再作傅里叶逆变换
?h??2?1U?x?h/2?=FFT-1?FFT?Uxexpi?*?????? ????2???2?非线性部分：
两边积分
?U?nU?pR(x)U ?zx?h1dU=??n?pR(x)?dz
xU?当h?0时
x?hxln最后有
U?x?h???n?pR?x??*h ??U?x?U?x?h??U?x?exp???n?pR?*h??
clear all
Lh=0; p=1;
omega=1; Dz=0;
%%------------------- n=2048; hx=0.06;
x=(-n/2:n/2-1)*hx; hw=2*pi/(n*hx);
w=fftshift((-n/2:n/2-1)*hw); %%------------------- q=exp(-1*(x).^2/2); % q=sech(x);
intensity=2;
u1(:,1)=(abs(q).^2)'; %-------------------
V=p*(cos(omega*x)).^2.*(1+Lh*exp(-x.^8/128)); %-------------------- L=500; nm=L*100; h=L/nm;
%------------------- for j=1:nm j;
D=exp(((i*w).^2/2)*h/2); qstep1=ifft(D.*fft(q));
n_index=ifft(fft(abs(qstep1).^2)./(Dz*w.^2+1)); N=exp((V+n_index)*h); qstep2=N.*qstep1; q=ifft(D.*fft(qstep2));
q=sqrt(intensity)*q/sqrt(sum(abs(q).^2)*hx); u=abs(q);
r=floor(2+(j-1)/L); u1(:,r)=u'; end
kin=-sum((q(3:end)-q(1:end-2)).^2)/4/hx;
p_i=sum(2*q.^2.*(abs(q).^2+V+n_index))*hx; b=(kin+p_i)/2/intensity z=0:L*h:L; figure(1)
mesh(x,z,u1'); view(0,90) figure(2)
plot(x,u1(:,end),'r',x,V,'b')
虚时间变换
?q1?2q2i??qq?pR?x?q?0 2?z2?x
做虚时间变换z??i?，则有：
?q?q1?q?q????i?z???i??i????
即
?q1?2q2???qq?pR?x?q?0 2??2?x
再将?记为z，有
?q1?2q2???qq?pR?x?q?0 2?z2?x
