分步傅里叶算法的MATLAB程序实现

分步傅里叶算法的MATLAB程序实现举例

模型：

线性部分： 非线性部分：2

i U1 U z

2 x

2 nU 0

2

2

n U x,z d z

n x

2其中：

d z exp z 或d z exp z

2

U z

=

i U2 x

2

两边同时对x变量作傅里叶变换

~

U z

=

i2

i

2

~

U

两边积分

x h21

~

U=x hix

~

U

x

2

i

2

dz

即

~

ln

U x h/2 ~

=

iU x

2

i

2

*

h2

最后有

~

~

U x h/2 =U x exp i

2h i * 2

2

再对x变量作作傅里叶逆变换

U x h/2 =FFT-1

FFT i2h U x exp 2 i *

2

U z

inU

两边积分

x h1 hx

U

=

xx

indz

当h 0时

ln

U x h U x

inh

最后有

U x h U x exp inh

2

折射率部分： n U x,

两边同时对x变量作傅里叶变换

n=FFT

~~

2

nd 2

x

2

d z * i n

2

~

2

FFT

n= 2

1 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 x

2

nU pR(x)U 0 n x

22

n x,z d

2

MATLAB程序实现：

线性部分：

U z

=1 U2 x

22

两边同时傅里叶变换

~

U z

=

12

i

2

~

U

两边积分

即

x h2x

1

~

~

U

U=

x hx

12

i

2

dz

ln

U x h/2 U x

~

~

=

12

i

2

*

h2

最后有

h 2 1

U 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 hx

=

x hx

n pR(x) dz

ln

U 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);

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')