分步傅里叶算法的MATLAB程序实现举例
模型：
线性部分： 非线性部分：
2i?U1?U?z?2?x2?nU?02?2
n?U?x,z??d?z?n?x2其中：
d?z??exp??z?或d?z??exp???z?
2
?U?z=i?U2?x2
两边同时对x变量作傅里叶变换
~?U?z=i2?i??2~U
两边积分
?x?h21~dU=x?h2ix~U?x2?i??2dz
即
~lnU?x?h/2?~=iU?x?2?i??2*h2
最后有
~~U?x?h/2?=U?x?exp?i2h???i??*?22?
?再对x变量作作傅里叶逆变换
U?x?h/2?=FFT-1??FFT??i2h???U?x???exp?2?i??*??2????
?U?z?inU
两边积分
?x?h1?hxUdU=?xxindz
当h?0时
lnU?x?h?U?x??inh
最后有
U?x?h??U?x?exp?inh?
2折射率部分： n?U?x,?z?两边同时对x变量作傅里叶变换
n=FFT?U?~~2?nd??z2
?x2??d?z?*?i??n ?2~2FFT?U???n= 21??d?z?再对x变量作作傅里叶逆变换
?FFT?U2?????-1?n=FFT?? 2?1??d?z????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
2得到
??U?z?1?U2?x2?nU?pR(x)U?0?n?x22
n?U?x,z??d2MATLAB程序实现：
线性部分：
?U?z=1?U2?x22
两边同时傅里叶变换
~?U?z=12?i??2~U
两边积分
?即
x?h2x1~~UdU=?x?h2x12?i??2dz
lnU?x?h/2?U?x?~~=12?i??2*h2
最后有
h?2?1U?x?h/2?=U?x?exp??i??*?
2??2~~再作傅里叶逆变换
h??2?1-1?U?x?h/2?=FFT?FFT?Uxexpi?*?????? ???2?2????非线性部分：
两边积分
?U?z1U?nU?pR(x)U
?当h?0时
x?hxdU=?x?hx?n?pR(x)?dz
lnU?x?h?U?x????n?pR?x???*h
最后有
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);
…… 此处隐藏0字 ……
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')