色散性的FDTD仿真求助

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [0lO0ik>G  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !UHX?<3r  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TPmb]j  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% _gQ_ixu  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4ULdf|oP"  
%%%%%%%%%初始化 95Q^7oI  
clear; cXK.^@du  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5YLc4z*  
%%%%%%%%%%%%%%%%系统参数 y<PQ$D)  
C)ic;!$Qhb  
TimeT=200;%迭代次数 V6_~"pRR=  
KE=2000;%网格树木 gSkY c{b  
kc=450;%源的位置 1V-sibE  
kpstart=500;%等离子体开始位置 e8{!Kjiz  
kpstop=1000;%等离子体终止位置 9>;CvR  
DispE=zeros(1,TimeT); vJe c+a  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% S7b7zJ8A  
%%%%%%%%%%%%%物理参数 _z>%h>L|g  
c0=3e8;%真空中波速 OV`li#H  
f=150e6; fzUG1|$e  
f1=75e6; Nb)Mh  
lamda=c0/f; XoGOY|2`6  
WL=50; 7,Y+FZ  
OMIGA=pi/WL; ptlag&Z  
zdelta=lamda/WL;%网格大小 fYlqaO4[  
dt=zdelta/(2*c0);%时间间隔 pW_mS|  
R,8Tt!n  
G)&!f)6  
u0=1e7;%碰撞频率 y 7z)lBy\  
fpe=1e7;%等离子体频率 ;WG
%)^e  
wpe=2*pi*fpe;%等离子体圆频率 (+(bw4V/  
epsz=1/(4*pi*9*10^9); % 真空介电常数 .dj}y jd]f  
mu=1/(c0^2*epsz);%磁常数 \zhCGDm1_  
ex_low_m1=0; &;U F,  
ex_low_m2=0; !{vZvy"  
ex_high_m1=0; Zi<(>@z2  
ex_high_m2=0; PNH>LT^  
?=b#H6vs  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a0's6C  
%%%%%%%%%%%%%初始化电磁场 -`5L;cxwk4  
Ex=zeros(1,KE); }>_  
Ex_Pre=zeros(1,KE); fJ+4H4K  
Hy=zeros(1,KE); zpg5
12\y  
Hy_Pre=zeros(1,KE); =.qPjp_Qd  
Dx=zeros(1,KE); hHgH'  
Dx_Pre=zeros(1,KE); ol}}c6  
Sx1=zeros(1,KE); =]o2{d  
Sx2=zeros(1,KE); A@e!~  
Sx3=zeros(1,KE); Bnw^W_  
Sx=zeros(1,KE); |Rz}bsrZ  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Lm+E?Ca  
%%%%%%%%%%%%%%%%%%开始计算 #JR$RH  
for T=1:TimeT NT{'BJ  
    %%%保存前一时间的电磁场 S$/SFB$)~W  
    Ex_Pre=Ex; 8F/zrPG  
    Hy_Pre=Hy; dw'P =8d  
    %%%%中间差分计算Dx n<ecVFft  
    for i=2:KE I(<Trn  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); ;A0ZcgF  
    end n(_wt##wE~  
    %%%%%%%%加入源 " "S&zN  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); le6eorK8  
       ;8w CQ  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); cbW=kQc_  
    if T<WL $?On,U  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 lU.aDmy<  
    end HBA|NV3.  
    %%%计算电场Ex JY%l1:}G3  
    for i=1:kpstart-1 3gv?rJV  
        Ex(i)=Dx(i)/epsz; eh,_g.  
    end AR c  
    for i=kpstop+1:KE w I[Hoi V  
        Ex(i)=Dx(i)/epsz; C|]Zpn#{K  
    end V!s#xXD}  
    if T>50 ?gCP"~  
        gw=0; a!.Y@o5Ku  
    end 
f'Rq#b@  
    i=kpstart; B[5<&  
         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); *r ('A  
         Ex(i)=Dx(i)/epsz-Sx1(i); df/7u}>9  
         Ex(i)=Ex(i); Rp>%umDyL  
     7wHd*{^9N  
    for i=kpstart+1:kpstop TJ|do`fw>  
         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); a..LbQQ  
         Ex(i)=Dx(i)/epsz-Sx1(i); h<7@3Ur  
         %Ex(i)=Ex(i); qR kPl!5  
    end \"d\b><R  
    Sx3=Sx2; 10_>EY`  
    Sx2=Sx1; sTvw@o*  
    Ex(1)=ex_low_m2; ,M :j5  
    ex_low_m2=ex_low_m1; p{&o{+c  
    ex_low_m1=Ex(2); 
'tt4"z2  
zL3I!& z2  
    Ex(KE)=ex_high_m2; A(z m  
    ex_high_m2=ex_high_m1; QiaB
ZAol  
    ex_high_m1=Ex(KE-1); sHQO*[[  
    %%%%%%%%%%%%%%%%%%计算磁场 `:d\L H  
    for i=1:KE-1
#;h> x  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); [Bj\h7G
 
    end vnqLcNB H  
    DispE(T)=Ex(kpstart+50); zek>]l`!  
    %if round(T/100)*100==T z}Vg4\x&  
      plot(Ex); hJ?PV@xy
 
      grid on; /1eeNbd  
      pause(0.00001); Ac\e>N  
      T #<es>~0!  
      % end 4W=fQx]  
     T%E/k# )q  
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ]hRs -x  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R_2#7Xs  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% h!tg+9%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% z?GtC{L9  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% iF_u/#  
%%%%%%%%%初始化 <
SdOb#2  
clear; bAUYJPRpy  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% b#n  
%%%%%%%%%%%%%%%%系统参数 #%{  
TimeT=3000;%迭代次数 w^:@g~  
KE=2000;%网格树木 BM5)SgK  
kc=450;%源的位置 %VE FruM  
kpstart=500;%等离子体开始位置 [8vqw(2Tm(  
kpstop=1000;%等离子体终止位置 fc4jbPp:M  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Z{2QDjAI;  
%%%%%%%%%%%%%物理参数 TQOJN  
c0=3e8;%真空中波速 b,47 EJ}  
zdelta=1e-9;%网格大小 Wxgs66  
dt=zdelta/(2*c0);%时间间隔 @KJmNM1]V  
f=900e12;%Gause脉冲的载频 w9Nk8OsL  
d=3e-15%脉冲底座宽度 aM:tg1g  
t0=2.25/f;%脉冲中心时间 !/1~  
u0=57e12%碰撞频率 [C!m,4  
fpe=2000e12;%等离子体频率 pV1;gqXNS  
wpe=2*pi*fpe;%等离子体圆频率 F*Ul#yX  
epsz=1/(4*pi*9*10^9); % 真空介电常数 Z=l2Po n  
mu=1/(c0^2*epsz);%磁常数 q?8#D  
ex_low_m1=0; '1d0 *5+6k  
ex_low_m2=0; 2o7o~r  
ex_high_m1=0; J ayax]u7J  
ex_high_m2=0; aHPx'R  
a0=2*u0/dt+(2/dt)^2; 2#8PM-3"  
a1=-8/(dt)^2; 4-9cp=\PE  
a2=-2*u0/dt+(2/dt)^2; (l/i#  
b0=wpe^2+2*u0/dt+(2/dt)^2; ;zSV~G6-  
b1=2*wpe^2-8/(dt)^2; ~/Y8wxg  
b2=wpe^2-2*u0/dt+(2/dt)^2; ClufP6'  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4$4Tx9C  
%%%%%%%%%%%%%初始化电磁场 zLPCWP.u  
Ex=zeros(1,KE); N9vNSmm  
Ex_Pre=zeros(1,KE); |T-Ytuy8  
Hy=zeros(1,KE); >b~Q%{1  
Hy_Pre=zeros(1,KE); W k"_lJ  
Dx=zeros(1,KE); $?9u;+jIR  
Dx_Pre=zeros(1,KE); rMJ4w['J=  
Sx1=zeros(1,KE);
w1
A-_  
Sx2=zeros(1,KE); V'9OGn2v  
Sx3=zeros(1,KE); 8jiBLZkRf  
Sx=zeros(1,KE); z_). -  
Dx=Ex; Djf~8q V!  
Dx1=Ex; I]~s{I(EK  
Dx2=Ex; 89W8c
J$yW  
Dx3=Ex; ^}kYJvqA  
Ex1=Ex; )}k"7"  
Ex2=Ex; QRF:6bAxsL  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% )P|/<>z  
%%%%%%%%%%%%%%%%%%开始计算 \2Kl]G(w%y  
for T=1:TimeT *48LQzc  
    %%%保存前一时间的电磁场 S)"vyGv  
    Ex_Pre=Ex; 3
Q,p,  
    Hy_Pre=Hy; bHzZ4i  
    %%%%中间差分计算Dx [7[$P.MS{  
    for i=2:KE 6K >(n  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); \\ZhM  
    end }}wSns  
    %%%%%%%%加入源 [p)2!]y  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); \P
*%u  
     ) urUaE  
    %%%计算电场Ex \;+b1  
    for i=1:kpstart-1 N8!e(YK_  
        Ex(i)=Dx(i)/epsz; YP{mzGdE&  
    end 9~ifST\  
    for i=kpstop+1:KE Rb <{o8  
        Ex(i)=Dx(i)/epsz; Q?'Ax"$D  
    end u@zT~\ h*  
    Dx3=Dx2; Mi`t$hmP  
Dx2=Dx1; - 
@KT#  
Dx1=Dx; rWBgYh  
    for i=kpstart:kpstop GF0Utp:Zf;  
                  Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); vVo# nzeZ5  
    end {s?hXB  
    Ex2=Ex1; Xg}~\|n  
Ex1=Ex; T n,Ifo3  
    Sx3=Sx2; W\e!rq  
    Sx2=Sx1; l_fERp#y  
    Ex(1)=ex_low_m2; Y1BxRd?D  
    ex_low_m2=ex_low_m1; RC!T1o~L  
    ex_low_m1=Ex(2); P[e#j  
L)Da1<O  
    Ex(KE)=ex_high_m2; iT1HbAT]  
    ex_high_m2=ex_high_m1; sZokiFJ  
    ex_high_m1=Ex(KE-1); 5jTA6s9zA  
    %%%%%%%%%%%%%%%%%%计算磁场 [VouG{  
    for i=1:KE-1 *C>
N  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); 3W3d $  
    end ! n?j)p.  
    plot(Ex); 1i Q(q\%  
    grid on; OW8TiM mK  
    pause(0.01); QFhQfn  
end