色散性的FDTD仿真求助

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% :^?-bppYW  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >wqWIw.w>  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >V)#y$Z  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ,|$1(z*a{c  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -:_3N2U=+  
%%%%%%%%%初始化 /!Ag/SmS!9  
clear; _X?_|!;J  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% '>|
Kd{J0  
%%%%%%%%%%%%%%%%系统参数 [AFR \{  
Mj>QV(L8t  
TimeT=200;%迭代次数 !U4YA1>>  
KE=2000;%网格树木 ECL{`m(#n  
kc=450;%源的位置 ]lGkZyUhI  
kpstart=500;%等离子体开始位置 z[[qrR  
kpstop=1000;%等离子体终止位置 8SroA$^n  
DispE=zeros(1,TimeT); (kX:@9Pn  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% :dipk,b?n  
%%%%%%%%%%%%%物理参数 u8?$W%eW  
c0=3e8;%真空中波速 6W YVHG  
f=150e6; <^> nR3E  
f1=75e6; 5#f_1 V  
lamda=c0/f; 3<)][
<Ud  
WL=50; Y]6dYq{k  
OMIGA=pi/WL; {Q$8p2W  
zdelta=lamda/WL;%网格大小 gA
EB  
dt=zdelta/(2*c0);%时间间隔 = =pQ V[  
\25EI]  
e`LvHU_0  
u0=1e7;%碰撞频率 =[aiW|Y  
fpe=1e7;%等离子体频率 E]vox~xK>  
wpe=2*pi*fpe;%等离子体圆频率 %?V~7tHm>  
epsz=1/(4*pi*9*10^9); % 真空介电常数 oc-&}R4=  
mu=1/(c0^2*epsz);%磁常数 +}udIi3:l  
ex_low_m1=0; voRb>xF  
ex_low_m2=0; /z!y[ri+J  
ex_high_m1=0; `j
'1V1  
ex_high_m2=0; 6X'0 T}  
%8DI)n#H  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EY !o#m  
%%%%%%%%%%%%%初始化电磁场 |ul{d|  
Ex=zeros(1,KE);
'tMD=MH  
Ex_Pre=zeros(1,KE); um/F:rp  
Hy=zeros(1,KE); Y#9bM$x7  
Hy_Pre=zeros(1,KE); Y<Ae_yLa  
Dx=zeros(1,KE); 3hJ51=_0^  
Dx_Pre=zeros(1,KE); yZ=O+H  
Sx1=zeros(1,KE); *7*cWO=  
Sx2=zeros(1,KE); {C]tS5$Z  
Sx3=zeros(1,KE); $rQ7"w J  
Sx=zeros(1,KE); 7ieAd/:_  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% H0B=X l[  
%%%%%%%%%%%%%%%%%%开始计算 Ed0}$b  
for T=1:TimeT 'HV@i)h0%V  
    %%%保存前一时间的电磁场 d7X7_  
    Ex_Pre=Ex; }vP
(SF6  
    Hy_Pre=Hy; \L?A4Qx)_  
    %%%%中间差分计算Dx =s.0 f:(  
    for i=2:KE '"Uhw$#t  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); )J
yB  
    end \O~/^ Y3U!  
    %%%%%%%%加入源 EUrIh2.Z  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2);
T%"wz3~  
       w?*79 u  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); |DsT $~D  
    if T<WL  Jc]k\U  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 ^(T~Qp  
    end /ioBc}]  
    %%%计算电场Ex _@)-#7  
    for i=1:kpstart-1 >*B59+1P  
        Ex(i)=Dx(i)/epsz; dqBN_P%  
    end yfqe6-8U  
    for i=kpstop+1:KE Fku<|1}&y  
        Ex(i)=Dx(i)/epsz; l%0-W  
    end @+Nf@LJ  
    if T>50 TntTR"6aD  
        gw=0; qM3NQ8Rm  
    end <7Yh<(R e^  
    i=kpstart; T_wh)B4xW  
         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); VS`
{k^^  
         Ex(i)=Dx(i)/epsz-Sx1(i); lOb(XH
9  
         Ex(i)=Ex(i); I%Z=O=  
     Hv' OO@z  
    for i=kpstart+1:kpstop 3TV4|&W;  
         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); @|w/`!}9q  
         Ex(i)=Dx(i)/epsz-Sx1(i); inq {" 6  
         %Ex(i)=Ex(i); gE*7[*2?t  
    end @Wm:Rz  
    Sx3=Sx2; l*CCnqE  
    Sx2=Sx1; Y HS/|-  
    Ex(1)=ex_low_m2; -#4QY70H t  
    ex_low_m2=ex_low_m1; }[c.OJ:   
    ex_low_m1=Ex(2); b,jo94.G  
|eye) E:  
    Ex(KE)=ex_high_m2; GKUjtPu  
    ex_high_m2=ex_high_m1; m#_M"B.cm  
    ex_high_m1=Ex(KE-1); G<rAM+B*g  
    %%%%%%%%%%%%%%%%%%计算磁场  (t@!0_5  
    for i=1:KE-1 fV6ddh  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); kbH@h2Ww  
    end g%ys|  
    DispE(T)=Ex(kpstart+50); BSr#;;\  
    %if round(T/100)*100==T hhU\$'0B-  
      plot(Ex); p,(W?.ZDN?  
      grid on; ;5T}@4m|r  
      pause(0.00001); vq3:N'  
      T #Rs5W  
      % end C(>g4.-p8  
     I;?PDhDb  
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% p`52  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% g1l:k1\Ht  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% G$CSZrP.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ,u-9e4  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ~u[1Vz4#3  
%%%%%%%%%初始化 mGmZ}H'{  
clear; "W9z>ezp  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% W,!7_nl"u  
%%%%%%%%%%%%%%%%系统参数 ~~t>;  
TimeT=3000;%迭代次数 f3.oc9G  
KE=2000;%网格树木 E2`9H-6e  
kc=450;%源的位置 !#e+!h@  
kpstart=500;%等离子体开始位置 :eevc7  
kpstop=1000;%等离子体终止位置 4((p?jbC  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Li 9$N"2  
%%%%%%%%%%%%%物理参数 A\E ))b9+  
c0=3e8;%真空中波速 &#.>-D{  
zdelta=1e-9;%网格大小 Z;0<k;#T(p  
dt=zdelta/(2*c0);%时间间隔 pDhUD}1G  
f=900e12;%Gause脉冲的载频 ?UeV5<TewS  
d=3e-15%脉冲底座宽度 l!\~T"-7;:  
t0=2.25/f;%脉冲中心时间 aJ1{9 5ea  
u0=57e12%碰撞频率 [I<J6=  
fpe=2000e12;%等离子体频率 8R(l~  
wpe=2*pi*fpe;%等离子体圆频率 yKm6 8n^  
epsz=1/(4*pi*9*10^9); % 真空介电常数 0|=y#`;,Z  
mu=1/(c0^2*epsz);%磁常数 7pnlS*E.  
ex_low_m1=0; @2_E9{T  
ex_low_m2=0; cw0uLMqr`  
ex_high_m1=0; K,*z8@  
ex_high_m2=0; REsw=P!b  
a0=2*u0/dt+(2/dt)^2; {QOy' 8/  
a1=-8/(dt)^2; P/q] u  
a2=-2*u0/dt+(2/dt)^2; g$/7km{TP  
b0=wpe^2+2*u0/dt+(2/dt)^2; ) <w`:wD  
b1=2*wpe^2-8/(dt)^2; EEx:Xk%5hX  
b2=wpe^2-2*u0/dt+(2/dt)^2; L[`8 :}M  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "zqa:D26  
%%%%%%%%%%%%%初始化电磁场 (_gt!i{h  
Ex=zeros(1,KE); cveQ6 -`K  
Ex_Pre=zeros(1,KE); ;!3: 3;  
Hy=zeros(1,KE); %):pfM;b  
Hy_Pre=zeros(1,KE); l=$?#^^ /  
Dx=zeros(1,KE); )%8st'  
Dx_Pre=zeros(1,KE); +4[9Eb'k=  
Sx1=zeros(1,KE); L E
FLKC  
Sx2=zeros(1,KE); |S:erYE,G  
Sx3=zeros(1,KE); 3RFU  
Sx=zeros(1,KE); +u&3pK>f  
Dx=Ex; "Rtt~["%  
Dx1=Ex; 6%wlz%Fp  
Dx2=Ex; =&FaMR2  
Dx3=Ex; W!+=`[Ff  
Ex1=Ex; lWP]}Uy=5~  
Ex2=Ex; @zLyG#kHY  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 'S&Zq:  
%%%%%%%%%%%%%%%%%%开始计算 hyhm{RC?[  
for T=1:TimeT (W[]}k;  
    %%%保存前一时间的电磁场 m6gMVon  
    Ex_Pre=Ex; Fp]ErDan  
    Hy_Pre=Hy; '0I>  
    %%%%中间差分计算Dx 'cc{sjG  
    for i=2:KE w2lO[o~x}  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); <;1M!.)5  
    end :C={Z}t/F  
    %%%%%%%%加入源 'UuHyC2Ha3  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); j|XL$Q  
     z+zEH9.'  
    %%%计算电场Ex q[+KQ,  
    for i=1:kpstart-1 HTm`_}G9  
        Ex(i)=Dx(i)/epsz; X7tBpyi  
    end 1v
Jj?Uqc  
    for i=kpstop+1:KE zpzxCzU  
        Ex(i)=Dx(i)/epsz; L{&Yh|}  
    end 2gEF$?+q?  
    Dx3=Dx2; h$%h w+"4  
Dx2=Dx1; WfTl\Dxw  
Dx1=Dx; ^D ;EbR  
    for i=kpstart:kpstop ?_T[]I'  
                  Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); )gAqWbkB  
    end WQBV~.<Yv  
    Ex2=Ex1;
&!wtH  
Ex1=Ex; 2dbn~j0  
    Sx3=Sx2; t7,$u-  
    Sx2=Sx1; g\iSc~%?  
    Ex(1)=ex_low_m2; /^X)>1)j  
    ex_low_m2=ex_low_m1; $IdU  
    ex_low_m1=Ex(2); >CHb;*U  
M<Dvhy[  
    Ex(KE)=ex_high_m2; t=jG$
A  
    ex_high_m2=ex_high_m1; # 00?]6`z  
    ex_high_m1=Ex(KE-1); 6k{gI.
SG  
    %%%%%%%%%%%%%%%%%%计算磁场 D^f;X.Qm  
    for i=1:KE-1 >Y|P+Z\7  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); ]A1'+!1$  
    end nXjSf  
    plot(Ex); e
]{=#  
    grid on; Ies` !W^  
    pause(0.01); )Z\Zw~L  
end

超过48小时，移动至相关板块继续讨论 =?N$0F!  
此帖锁定 +"|TPKas  
鼠标轻点这里即可参与讨论