色散性的FDTD仿真求助

小弟初学FDTD，最近在做一个光色散性的仿真，仿真的原始微分方程是简化的薛定谔方程，这个方程左边是关于z 的一阶微分右边是关于t 的二阶微粉，这样化成差分方程时就会出现∆z/∆ t^2 的项，若∆ t=∆z/c ，则∆z/∆ t^2 =c^2/∆z，这一项就受到仿真步长∆z很大的影响，有时会出现不收敛的情况，请问这种情况下的∆z应该如何选取?

薛定谔方程和MAXWELL方程有一些差异。薛定谔方程的差分离散方法不再是FDTD，FDTD是特指MAXWELL方程的差分离散格式。
对你的情况，“方程左边是关于z 的一阶微分右边是关于t 的二阶微分”，时间和空间的微分阶数不一样，属于PARABOLIC方程，∆ t=∆z/c是针对FDTD，对你的情况不适宜。对于你的现在的情况，在差分离散格式中要用∆z/∆ t^2 =c。
P.S., 薛定谔方程一般是“时间一阶微分”和“空间二阶微分”，你的方程好像不是正常的薛定谔方程。无论怎样，如果你的描述是正确的，上面描述的方法是适用于你的。

你的薛定谔方程怎么右边有二阶时间项，你可以把此方程的图片贴出来。
未知函数的空间、时间离散是独立的，不相关的。空间、时间步长的关系应由方程的稳定条件来决定，不同方程形式稳定条件不一样。

可以用FDTD了，其实就是个色散介质的问题，因为它的那个参数比如介电常数和频率有关，所以才多一相对时间的导数。
求解方法也有很多种，移位（位移?）算法，Z变化等等都可以的。
可惜我自己的电脑里面都是游戏，连office和pdf都没装，看不了以前的资料，要不然倒是可以上传上来。

一般在做差分时，f ' = f(n+1) - f(n) / delta_t; 实际上这是一个一阶近似。如果用exp来表示参数就可以解决这个问题了。

能说得再详细一点吗

这是冷等离子体的FDTD Z变换公式推导。你的应该更简单一些，毕竟系数中只有一次导数。

或者按照楼上的说法，你把公式贴出来，大家看看。

就是这个公式

我怎么看不见公式呢。你用图片传上来看看。或者弄成附件让我下一下。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%初始化
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%系统参数
TimeT=200;%迭代次数
KE=2000;%网格树木
kc=450;%源的位置
kpstart=500;%等离子体开始位置
kpstop=1000;%等离子体终止位置
DispE=zeros(1,TimeT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%物理参数
c0=3e8;%真空中波速
f=150e6;
f1=75e6;
lamda=c0/f;
WL=50;
OMIGA=pi/WL;
zdelta=lamda/WL;%网格大小
dt=zdelta/(2*c0);%时间间隔
u0=1e7;%碰撞频率
fpe=1e7;%等离子体频率
wpe=2*pi*fpe;%等离子体圆频率
epsz=1/(4*pi*9*10^9); % 真空介电常数
mu=1/(c0^2*epsz);%磁常数
ex_low_m1=0;
ex_low_m2=0;
ex_high_m1=0;
ex_high_m2=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%初始化电磁场
Ex=zeros(1,KE);
Ex_Pre=zeros(1,KE);
Hy=zeros(1,KE);
Hy_Pre=zeros(1,KE);
Dx=zeros(1,KE);
Dx_Pre=zeros(1,KE);
Sx1=zeros(1,KE);
Sx2=zeros(1,KE);
Sx3=zeros(1,KE);
Sx=zeros(1,KE);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%开始计算
for T=1:TimeT
    %%%保存前一时间的电磁场
    Ex_Pre=Ex;
    Hy_Pre=Hy;
    %%%%中间差分计算Dx
    for i=2:KE
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));
    end
    %%%%%%%%加入源
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2);
      
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T);
    if T<WL
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数
    end
    %%%计算电场Ex
    for i=1:kpstart-1
        Ex(i)=Dx(i)/epsz;
    end
    for i=kpstop+1:KE
        Ex(i)=Dx(i)/epsz;
    end
    if T>50
        gw=0;
    end
    i=kpstart;
         Sx1(i)=(1+exp(-1*u0*dt))*Sx2(i)-exp(-1*u0*dt)*Sx3(i)+(wpe^2*dt/(2*u0))*(1-exp(-1*u0*dt))*Ex(i);
         Ex(i)=Dx(i)/epsz-Sx1(i);
         Ex(i)=Ex(i);
    
    for i=kpstart+1:kpstop
         Sx1(i)=(1+exp(-1*u0*dt))*Sx2(i)-exp(-1*u0*dt)*Sx3(i)+(wpe^2*dt/(u0))*(1-exp(-1*u0*dt))*Ex(i);
         Ex(i)=Dx(i)/epsz-Sx1(i);
         %Ex(i)=Ex(i);
    end
    Sx3=Sx2;
    Sx2=Sx1;
    Ex(1)=ex_low_m2;
    ex_low_m2=ex_low_m1;
    ex_low_m1=Ex(2);
    Ex(KE)=ex_high_m2;
    ex_high_m2=ex_high_m1;
    ex_high_m1=Ex(KE-1);
    %%%%%%%%%%%%%%%%%%计算磁场
    for i=1:KE-1
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i));
    end
    DispE(T)=Ex(kpstart+50);
    %if round(T/100)*100==T
      plot(Ex);
      grid on;
      pause(0.00001);
      T
      % end
    
end

这是一维等离子体FDTD源代码，Matlab程序，其中的介电系数为epsonr= epson0 * ( 1 - wpe^2 / (w * ( w - i * up))) 。所以你会发现在求解的时候用了exp( - u0 * dt) ，而不是1 - u0 * dt ，这就是我前面说的，其实1 - u0 * dt 是exp(-u0*dt)的一阶近似。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%初始化
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%系统参数
TimeT=3000;%迭代次数
KE=2000;%网格树木
kc=450;%源的位置
kpstart=500;%等离子体开始位置
kpstop=1000;%等离子体终止位置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%物理参数
c0=3e8;%真空中波速
zdelta=1e-9;%网格大小
dt=zdelta/(2*c0);%时间间隔
f=900e12;%Gause脉冲的载频
d=3e-15%脉冲底座宽度
t0=2.25/f;%脉冲中心时间
u0=57e12%碰撞频率
fpe=2000e12;%等离子体频率
wpe=2*pi*fpe;%等离子体圆频率
epsz=1/(4*pi*9*10^9); % 真空介电常数
mu=1/(c0^2*epsz);%磁常数
ex_low_m1=0;
ex_low_m2=0;
ex_high_m1=0;
ex_high_m2=0;
a0=2*u0/dt+(2/dt)^2;
a1=-8/(dt)^2;
a2=-2*u0/dt+(2/dt)^2;
b0=wpe^2+2*u0/dt+(2/dt)^2;
b1=2*wpe^2-8/(dt)^2;
b2=wpe^2-2*u0/dt+(2/dt)^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%初始化电磁场
Ex=zeros(1,KE);
Ex_Pre=zeros(1,KE);
Hy=zeros(1,KE);
Hy_Pre=zeros(1,KE);
Dx=zeros(1,KE);
Dx_Pre=zeros(1,KE);
Sx1=zeros(1,KE);
Sx2=zeros(1,KE);
Sx3=zeros(1,KE);
Sx=zeros(1,KE);
Dx=Ex;
Dx1=Ex;
Dx2=Ex;
Dx3=Ex;
Ex1=Ex;
Ex2=Ex;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%开始计算
for T=1:TimeT
    %%%保存前一时间的电磁场
    Ex_Pre=Ex;
    Hy_Pre=Hy;
    %%%%中间差分计算Dx
    for i=2:KE
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));
    end
    %%%%%%%%加入源
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2);
    
    %%%计算电场Ex
    for i=1:kpstart-1
        Ex(i)=Dx(i)/epsz;
    end
    for i=kpstop+1:KE
        Ex(i)=Dx(i)/epsz;
    end
    Dx3=Dx2;
Dx2=Dx1;
Dx1=Dx;
    for i=kpstart:kpstop
                  Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i));
    end
    Ex2=Ex1;
Ex1=Ex;
    Sx3=Sx2;
    Sx2=Sx1;
    Ex(1)=ex_low_m2;
    ex_low_m2=ex_low_m1;
    ex_low_m1=Ex(2);
    Ex(KE)=ex_high_m2;
    ex_high_m2=ex_high_m1;
    ex_high_m1=Ex(KE-1);
    %%%%%%%%%%%%%%%%%%计算磁场
    for i=1:KE-1
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i));
    end
    plot(Ex);
    grid on;
    pause(0.01);
end