分步傅里叶算法的MATLAB程序实现
一、模型:
其中:
或
线性部分:
两边同时对x变量作傅里叶变换
两边积分
即
最后有
再对x变量作作傅里叶逆变换
非线性部分:
两边积分
当h->0时
最后有
折射率部分:
两边同时对x变量作傅里叶变换
再对x变量作作傅里叶逆变换
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')
二、虚时间变换
作虚时间变换:
得到
算法实现:
线性部分:
两边同时傅里叶变换
两边积分
即
最后有
再作傅里叶逆变换
非线性部分:
两边积分
当h->0时
最后有
程序源码如下:
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')
虚时间变换
做虚时间变换
则有:
即
再将记为,有