色散性的FDTD仿真求助

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (R Ttz  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% <.
RgMPi  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% :?yv0Iu  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0*,r  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5C2 *f4|  
%%%%%%%%%初始化 ?%;)> :3N  
clear; 8n1Sy7K!;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HR)joD*q;[  
%%%%%%%%%%%%%%%%系统参数 Rs5G5W@"A  
!QP~#a%  
TimeT=200;%迭代次数 !.x(lOqf  
KE=2000;%网格树木 lvk(q\-f  
kc=450;%源的位置 fWF\V[  
kpstart=500;%等离子体开始位置 EPdR-dC^wE  
kpstop=1000;%等离子体终止位置 wcB-)Ra  
DispE=zeros(1,TimeT); +9,"ne1'e  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ~@TNVkw  
%%%%%%%%%%%%%物理参数 m9xO& @#vx  
c0=3e8;%真空中波速 v0\2%PC  
f=150e6; xZbm,.v  
f1=75e6; ZZ?=^g  
lamda=c0/f; t\u0\l>  
WL=50; ;l+3l ez  
OMIGA=pi/WL; s+o/:rrxY  
zdelta=lamda/WL;%网格大小 [
@&  
dt=zdelta/(2*c0);%时间间隔 tAc[r)xFw  
Sxf<8Px9i  
d];E99}  
u0=1e7;%碰撞频率 :] {+
3A  
fpe=1e7;%等离子体频率 gb8nST$r  
wpe=2*pi*fpe;%等离子体圆频率 VO[s:e9L  
epsz=1/(4*pi*9*10^9); % 真空介电常数 uu]<R@!J  
mu=1/(c0^2*epsz);%磁常数 FZb\VUmnV  
ex_low_m1=0; IP/ zFbc  
ex_low_m2=0; w6v P a  
ex_high_m1=0; RcMW%q$dG  
ex_high_m2=0; M|6 W<y  
ep>S$a*|  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |jF)~k6  
%%%%%%%%%%%%%初始化电磁场 F;5.nKo  
Ex=zeros(1,KE); ,W7\AY07]  
Ex_Pre=zeros(1,KE); 4LJUO5(y@  
Hy=zeros(1,KE); l 0jjLqm:  
Hy_Pre=zeros(1,KE); DRj\i6-v  
Dx=zeros(1,KE); )z$VQ=]"  
Dx_Pre=zeros(1,KE); 4Poi:0oOys  
Sx1=zeros(1,KE); * t{A=Wk  
Sx2=zeros(1,KE); `A$yF38!  
Sx3=zeros(1,KE); kbBX\*{yh  
Sx=zeros(1,KE); \0'o*nlJ  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% L]9uY
  
%%%%%%%%%%%%%%%%%%开始计算 ko$bCG%  
for T=1:TimeT M($dh9A_  
    %%%保存前一时间的电磁场 ,8cw jS2E  
    Ex_Pre=Ex; ;v*$6DIC5  
    Hy_Pre=Hy; bu,Z'  
    %%%%中间差分计算Dx Vv0dBFe  
    for i=2:KE &ZClv"6  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); |:L<Ko  
    end <s)+V6\E  
    %%%%%%%%加入源 &tkkn2t  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); ]?/[& PP,  
       4:y;<8+j\  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); W9t%:wF  
    if T<WL sq*d?<:3  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 JyvXNV,  
    end QjVP]C}p  
    %%%计算电场Ex =YR/X@&  
    for i=1:kpstart-1 ''s]6Jjw  
        Ex(i)=Dx(i)/epsz; GGo n
A  
    end =Bu d!  
    for i=kpstop+1:KE \\Fl,'  
        Ex(i)=Dx(i)/epsz; c{ +Y$  
    end [^
A93F  
    if T>50 :%{MMhbx  
        gw=0; pMHY2t  
    end 6,ylkf3  
    i=kpstart; %19TJn%J$  
         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); TWK(vEDM  
         Ex(i)=Dx(i)/epsz-Sx1(i); 8Tyf#`'I  
         Ex(i)=Ex(i); bT@3fuL4  
     EXK~Zf|&Z  
    for i=kpstart+1:kpstop Ha)e
eE$  
         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); (B.J8`h }  
         Ex(i)=Dx(i)/epsz-Sx1(i); q UY;CEf  
         %Ex(i)=Ex(i); V)^nVD)e  
    end oQBfDD0  
    Sx3=Sx2; wwF]+w%lOw  
    Sx2=Sx1; q%]0%S?  
    Ex(1)=ex_low_m2; _1VtVfiZ{  
    ex_low_m2=ex_low_m1; D[x0sly  
    ex_low_m1=Ex(2); 7N+No.vR.  
H^Ik FEVs  
    Ex(KE)=ex_high_m2; 
.36z  
    ex_high_m2=ex_high_m1; 
g|a2z_R  
    ex_high_m1=Ex(KE-1); <VS\z(K  
    %%%%%%%%%%%%%%%%%%计算磁场 MW=2GhD=  
    for i=1:KE-1 Q|rrbxb  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); ?{mFQ
  
    end .Vj;[p8  
    DispE(T)=Ex(kpstart+50); -5l74f!i  
    %if round(T/100)*100==T U~azI(1"W  
      plot(Ex); wLAGe'GX  
      grid on; {n2mh%I  
      pause(0.00001); ,9^wKS!7$  
      T we9R4*j  
      % end [RiCa  
     Z)Nl\e& M  
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >IT19(J;A  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (=\))t8J  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *#y9Pve  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% `<#Ufi*c  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% A )q=.C#e  
%%%%%%%%%初始化 qpEK36Js  
clear; z JBcz,  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uUIjntSF(  
%%%%%%%%%%%%%%%%系统参数 ._X|Ye9/  
TimeT=3000;%迭代次数 !_P-?u  
KE=2000;%网格树木 ,
tEd>  
kc=450;%源的位置 jtH>&O  
kpstart=500;%等离子体开始位置 .EfGL_  
kpstop=1000;%等离子体终止位置  8MZ:=  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (ah^</  
%%%%%%%%%%%%%物理参数 }+/F?_I= %  
c0=3e8;%真空中波速 "AuU5G 9'I  
zdelta=1e-9;%网格大小 W Te1E,M  
dt=zdelta/(2*c0);%时间间隔 oF(=@UL  
f=900e12;%Gause脉冲的载频 yDORL| E'  
d=3e-15%脉冲底座宽度 hY(q@_s  
t0=2.25/f;%脉冲中心时间 SHA6;y+U/~  
u0=57e12%碰撞频率 [G<SAWFg7  
fpe=2000e12;%等离子体频率 j>I.d+  
wpe=2*pi*fpe;%等离子体圆频率 3vc2t6S%*  
epsz=1/(4*pi*9*10^9); % 真空介电常数 8ioxb`U  
mu=1/(c0^2*epsz);%磁常数 ETQL,t9m  
ex_low_m1=0; .L=C7w1  
ex_low_m2=0; zI&).  
ex_high_m1=0; _8{6&AmIw  
ex_high_m2=0; m\"X%Y#  
a0=2*u0/dt+(2/dt)^2; gyT3[*eh  
a1=-8/(dt)^2; W*Gp0pX  
a2=-2*u0/dt+(2/dt)^2; >U#j\2!Sg  
b0=wpe^2+2*u0/dt+(2/dt)^2; oFDJwOJ'Bj  
b1=2*wpe^2-8/(dt)^2; eFz!`a^dX  
b2=wpe^2-2*u0/dt+(2/dt)^2; {SJnPr3R  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 >:RFCo  
%%%%%%%%%%%%%初始化电磁场 BnPL>11Y  
Ex=zeros(1,KE); o#frNT}  
Ex_Pre=zeros(1,KE); ,0^9VWZV  
Hy=zeros(1,KE); w<me(!-'  
Hy_Pre=zeros(1,KE); Lv<)Dur0K  
Dx=zeros(1,KE); 3IYbgUG  
Dx_Pre=zeros(1,KE); 2O+
fjs  
Sx1=zeros(1,KE); "*oN~&flc  
Sx2=zeros(1,KE); 
`XK+Y  
Sx3=zeros(1,KE); ^!x}e+ o  
Sx=zeros(1,KE); 1vL$k[^&d  
Dx=Ex; NVG`XL  
Dx1=Ex; &aWY{ ?_  
Dx2=Ex; *OR(8;  
Dx3=Ex; v>'mW  
Ex1=Ex; 1g1gu=|Q  
Ex2=Ex; U\`yLsKvH`  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% F94Qb}  
%%%%%%%%%%%%%%%%%%开始计算 9 Xx4,#?  
for T=1:TimeT kOLS<>.  
    %%%保存前一时间的电磁场 idGhWV'  
    Ex_Pre=Ex; H
\RuYCn2G  
    Hy_Pre=Hy; fudLm  
    %%%%中间差分计算Dx gt:Ot0\7  
    for i=2:KE Xb5$ijH  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); bl-t>aO*.V  
    end 1`@rAA>h'  
    %%%%%%%%加入源 v}^ f8nVR  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); r/BiR0$E  
     ^ ^R4%C  
    %%%计算电场Ex c69M   
    for i=1:kpstart-1 ko<VB#pOMr  
        Ex(i)=Dx(i)/epsz; ^`Qh*:T$  
    end liG3   
    for i=kpstop+1:KE a&~]77)  
        Ex(i)=Dx(i)/epsz; dD}!E  
    end Bl8&g]dk  
    Dx3=Dx2; "Qxn}$6-  
Dx2=Dx1; aFrVP  
Dx1=Dx; j,CMcP7A -  
    for i=kpstart:kpstop gKay3}w  
                  Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); ||vQW\g  
    end js8GK  
    Ex2=Ex1; $[-{Mm  
Ex1=Ex; *
3We5  
    Sx3=Sx2; 4,g3 c  
    Sx2=Sx1; S&m5]h!D  
    Ex(1)=ex_low_m2; X>6VucH{\  
    ex_low_m2=ex_low_m1;
OJ\rT.{  
    ex_low_m1=Ex(2); M[ZuXH}  
)B'U_*  
    Ex(KE)=ex_high_m2; Lu?)Rya  
    ex_high_m2=ex_high_m1; *tZ#^YG{(  
    ex_high_m1=Ex(KE-1); Q|W!m0XO  
    %%%%%%%%%%%%%%%%%%计算磁场 ,*$/2nB^  
    for i=1:KE-1 S.Fip_  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); 5`3f"(ay/  
    end ]Zf@NY  
    plot(Ex); Eh)VU_D  
    grid on; 7" wn024  
    pause(0.01); x{|n>3l`b9  
end