24小时热门版块排行榜

返回列表

zhang-zhen

金虫 (职业作家)

应助: 16 (小学生)
金币: 1191.1
散金: 2659
红花: 29
帖子: 3348
在线: 616.7小时
虫号: 1286891
注册: 2011-05-05
性别: GG
专业: 光学

[交流] 分布傅里叶算法求解非线性薛定谔方程报错

%有没有懂分布傅里叶算法的高手，帮忙检查下下面的非线性薛定谔方程求解代码究竟错在哪里？内循环里面的其它变量和jj的关系没有看到，不知道错误是否是出现在这里。谢谢：
%：）的符号用冒号和右括号代替即可
源代码如下：
clear all;clc;clf;
po=200;%输入光强，单位w
alpha=1.2;  %光纤损耗值，单位为db/km
gamma=20; %光纤非线性参数
to=1e-12;    %初始脉冲宽度，单位为秒
c=2;    %第一次计算输入的啁啾参数
% b2=1e-12*1e-12/0.1e3; %载波的倒数
b2=1000; %载波的倒数
cputime=0;
tic;
ln=1;
i=sqrt(-1);
pi=3.1415926535;
alph=alpha/(4.343);
belta=10;          %belta为色散系数，单位为ps/km,待定量
ld=(to^2)/(abs(b2)); %扩散长度，单位是m
ao=sqrt(po);          %光振幅
tau=-4096e-12:1e-12:4095e-12;
dt=1e-12;
h=1000;
for ii=0.1:0.1:1.5 %不同的光纤长度不同，这个量可变
z=ii*ld;
u=ao*exp(-(1+i*(-c))/2*(tau/to).^2);
figure(1)
plot(abs(u),'-r');
xlim([4090 4105]);
title('input pulse'); xlabel('time'); ylabel('amplitude');
% grid on;
hold on;
l=max(size(u));

fwhm1=find(abs(u)> abs(max(u)/2));
fwhm1=length(fwhm1);
dw=1/l/dt*2*pi;
w=(-1*l/2:1:l/2-l)*dw;
u=fftshift(u); %零延迟对中的谱
w=fftshift(w); %零延迟对中的谱
spectrum=fft(fftshift(u));%快速离散傅立叶变换
for jj=h:h:15*h
g1=(i*belta*w.^2-2*alph)/4*z;
spectrum=spectrum.*exp(g1) ; %g1为线性算符e的指数表达式
f=ifft(spectrum); %快速离散反傅立叶变换
g2=(i/2*belta*w.^2-alph)/2*z;
f=f.*exp(g2);%g2为非线性算符e的指数表达式
spectrum=fft(f); %快速离散傅立叶变换
spectrum=spectrum.*exp(g1) ;
end
f=ifft(spectrum); %快速离散反傅立叶变换
op_pulse(ln,

=abs(f);%保存在所有间隔点上的输出脉冲
fwhm=find(abs(f)>abs(max(f)/2));
fwhm=length(fwhm);
ratio=fwhm/fwhm1;
pbratio(ln)=ratio;
dd=atand((abs(imag(f)))/(abs(real(f))));
phadisp(ln)=dd;%保存脉冲相位
ln=ln+1;
end
toc;
cputime=toc;
figure(2)
mesh(op_pulse(1:1:ln-1,

);
title('pulse evolution');
xlim([4090 4105]);
xlabel('time'); ylabel('distance'); zlabel('amplitude');
figure(3)
plot(pbratio(1:1:ln-1),'k');
xlabel('number of steps');
ylabel('pulse broadening ratio');
grid on;
hold on;
figure(4)
plot(phadisp(1,1:1:ln-1),'k');
xlabel('distance travelled');
ylabel('phase change');
grid on;
hold on;
figure(5)
plot(f(1,

,'k');
xlabel('distance travelled');
ylabel('f');
disp('cpu time:'), disp(cputime)

回复此楼