[资料共享]一维等离子体FDTD的Matlab源代码(两种方法）

[post]一维等离子体FDTD的Matlab源代码 6 4da~SEn  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ^KJIT3J(#  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ZrFC#wJb  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% }]H_|V*f  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \ oIVE+L/P  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% X|7Y|0o  
%%%%%%%%%初始化 XK>/i}y  
clear; mSzBNvci  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -)tu$W*  
%%%%%%%%%%%%%%%%系统参数 #`mo5  
V4OhdcW{  
TimeT=200;%迭代次数 !s]LWCX+|  
KE=2000;%网格树木 [$Ld>`3  
kc=450;%源的位置 ~sQN\]5VW  
kpstart=500;%等离子体开始位置 /0mbG!Ac  
kpstop=1000;%等离子体终止位置 NVMhbpX6  
DispE=zeros(1,TimeT); h+x"?^  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ?1('s0s\,  
%%%%%%%%%%%%%物理参数 ~Cj55S+  
c0=3e8;%真空中波速 +M6qbIO  
f=150e6;  (Ia}]q  
f1=75e6; 451r!U1Z  
lamda=c0/f; hb"t8_--c  
WL=50; sgo({zA`i  
OMIGA=pi/WL; {7)D/WY5  
zdelta=lamda/WL;%网格大小 Q+[e)YO)  
dt=zdelta/(2*c0);%时间间隔 4cql?W(D  
cRX0i;zag  
^Q]*CU+C  
u0=1e7;%碰撞频率 ]}cai1  
fpe=1e7;%等离子体频率 fi%u]  
wpe=2*pi*fpe;%等离子体圆频率 j3rBEQ,R  
epsz=1/(4*pi*9*10^9); % 真空介电常数 $LZf&q:\]*  
mu=1/(c0^2*epsz);%磁常数 ]+W+8)f1M  
ex_low_m1=0; Nf>1`eP  
ex_low_m2=0; `av8|;  
ex_high_m1=0; l'(Cxhf.W  
ex_high_m2=0; *lg1iP{]  
a2*WZc
`  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @M?N[LG  
%%%%%%%%%%%%%初始化电磁场 3C8'0DB  
Ex=zeros(1,KE); 
Pn
5@7~  
Ex_Pre=zeros(1,KE); a4X J0Tm  
Hy=zeros(1,KE); dfe 9)m>  
Hy_Pre=zeros(1,KE); /
b20!3  
Dx=zeros(1,KE);  c/I.`@  
Dx_Pre=zeros(1,KE); Roy0?6O  
Sx1=zeros(1,KE); @YP\!#"8  
Sx2=zeros(1,KE); [SgP1>M  
Sx3=zeros(1,KE); ]?xF'3#  
Sx=zeros(1,KE); k'wF+>  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E)>~0jv  
%%%%%%%%%%%%%%%%%%开始计算 UnZ*"%  
for T=1:TimeT =eSG7QfS  
    %%%保存前一时间的电磁场 "TKf"z
c  
    Ex_Pre=Ex; [syuoJ  
    Hy_Pre=Hy; %L{H_;z  
    %%%%中间差分计算Dx dZRz'd  
    for i=2:KE `pN"T?Pk  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); l2=.;7IV  
    end X",fp  
    %%%%%%%%加入源 nbw&+dcJ8  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); yyrCO"eh  
        "=H7p3  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); Fm{Ri=X<:  
    if T<WL tsU.c"^n  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 GdR>S('  
    end ,v$gQU2  
    %%%计算电场Ex N- ?U2V  
    for i=1:kpstart-1 ,>2ijk#  
        Ex(i)=Dx(i)/epsz; J& +s  
    end t qbS!r  
    for i=kpstop+1:KE 1#Dpj.cO#  
        Ex(i)=Dx(i)/epsz; z['>`Kt  
    end H]Q Z4(  
    if T>50 $.cNY+ k  
        gw=0; c}Y(Myd  
    end Q}W6?XDu  
    i=kpstart; XY1NTo.=  
         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); O:RPH{D  
         Ex(i)=Dx(i)/epsz-Sx1(i); ,y3o ,gl  
         Ex(i)=Ex(i); J:'c
j5@  
     WO)rJr!C  
    for i=kpstart+1:kpstop 6t TLyI$+  
         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); ?y'KX]/  
         Ex(i)=Dx(i)/epsz-Sx1(i); KB7CO:  
         %Ex(i)=Ex(i); CY0|.x  
    end p(%7|'  
    Sx3=Sx2; Dz]&|5'N  
    Sx2=Sx1; DL|,:2`  
    Ex(1)=ex_low_m2; Tm_AoZH  
    ex_low_m2=ex_low_m1; sZPPS&KoP3  
    ex_low_m1=Ex(2); v
X)JJ|g  
R(=Lhz6R4  
    Ex(KE)=ex_high_m2; Vur$t^zE  
    ex_high_m2=ex_high_m1; &iR>:=ksN  
    ex_high_m1=Ex(KE-1); ly}6zOC\  
    %%%%%%%%%%%%%%%%%%计算磁场 yd`xm
c)  
    for i=1:KE-1 < .. ['sj'3cW-  
(X,Ua+{  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  N-`Vb0;N  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8@]*X,umc  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% D9,609w  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% $3<,"&;Ecs  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c[\ :^w^I6  
%%%%%%%%%初始化 F-[zuYGp  
clear; PtCO';9[  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T6m#sVq  
%%%%%%%%%%%%%%%%系统参数 ^!^6 |[  
TimeT=3000;%迭代次数 W2/FGJD  
KE=2000;%网格树木 (MhC83|?  
kc=450;%源的位置 F1)B-wW  
kpstart=500;%等离子体开始位置 2_M+akqy^  
kpstop=1000;%等离子体终止位置 ph{p[QI:{X  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HM57b>
6  
%%%%%%%%%%%%%物理参数 fZ1v|  
c0=3e8;%真空中波速 QA>(}u\+  
zdelta=1e-9;%网格大小 (XA=d 4  
dt=zdelta/(2*c0);%时间间隔 IhnHNY]<g  
f=900e12;%Gause脉冲的载频 b~fX=!M  
d=3e-15%脉冲底座宽度 Lh3>xZy"-z  
t0=2.25/f;%脉冲中心时间 ^ CVhV  
u0=57e12%碰撞频率 ?T=]?[  
fpe=2000e12;%等离子体频率 Wt5x*p-!C  
wpe=2*pi*fpe;%等离子体圆频率 Py7!_T
X  
epsz=1/(4*pi*9*10^9); % 真空介电常数 .w2QiJ  
mu=1/(c0^2*epsz);%磁常数 b?9c\-}  
ex_low_m1=0; #=F"PhiX`  
ex_low_m2=0; !`=ms1%U  
ex_high_m1=0; kms&o=^  
ex_high_m2=0; MJNY#v3  
a0=2*u0/dt+(2/dt)^2; Vcn04j#Q  
a1=-8/(dt)^2; u>c\J|K_V  
a2=-2*u0/dt+(2/dt)^2; ?~~sOf
AP  
b0=wpe^2+2*u0/dt+(2/dt)^2; Z0&^U#]  
b1=2*wpe^2-8/(dt)^2; 0q'd }DW  
b2=wpe^2-2*u0/dt+(2/dt)^2; >dKK [E/[d  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% b~DtaGh  
%%%%%%%%%%%%%初始化电磁场 9ZvBsG)  
Ex=zeros(1,KE); fm$eJu  
Ex_Pre=zeros(1,KE); vzV,} S*c  
Hy=zeros(1,KE); n][/c_]q  
Hy_Pre=zeros(1,KE); 3ThBy'
 
Dx=zeros(1,KE); j.FA
!4L  
Dx_Pre=zeros(1,KE); 4w,=6|#  
Sx1=zeros(1,KE); x,$N!X  
Sx2=zeros(1,KE); J-*&&  
Sx3=zeros(1,KE); 1Vq]4_09g1  
Sx=zeros(1,KE); lOIBX@K E  
Dx=Ex; mr:;Wwd  
Dx1=Ex; /2}o:vLj  
Dx2=Ex; Q#C;4)e  
Dx3=Ex; MuNM)pyxp  
Ex1=Ex; ]=\Mf<  
Ex2=Ex; }OY]mAv-B  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ' >(])Oq,  
%%%%%%%%%%%%%%%%%%开始计算 Z$qFjWp  
for T=1:TimeT 3TUW+#[Gu  
    %%%保存前一时间的电磁场 *Q2;bmIc  
    Ex_Pre=Ex; tHNvb\MR$  
    Hy_Pre=Hy; <$\vL  
    %%%%中间差分计算Dx n]M1'yU  
    for i=2:KE z_%G{H+:l  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); V3;4,^=6Dd  
    end `$og]Dn;  
    %%%%%%%%加入源 $=dp)  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); RjS;Ck@;  
     H /Idc,*  
    %%%计算电场Ex fDc>E+,  
    for i=1:kpstart-1 J ytY6HF  
        Ex(i)=Dx(i)/epsz; Rz}?@zh_8  
    end xdWfrm$;ZA  
    for i=kpstop+1:KE `5 py6,  
        Ex(i)=Dx(i)/epsz; w0Q
N5?  
    end 9hAS#|vK  
    Dx3=Dx2; [6x-c;H_4  
Dx2=Dx1; 0_yE74i  
Dx1=Dx; xe^*\6Y  
    for i=kpstart:kpstop gfQ&U@N  
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); +EJwWDJ!%  
    end #PnuR2s7.  
    Ex2=Ex1; `|K,E  
Ex1=Ex; =! v.VF\;  
    Sx3=Sx2; )6|7L)Dk  
    Sx2=Sx1; QS2J271E}  
    Ex(1)=ex_low_m2; Zu(eYH=Q  
    ex_low_m2=ex_low_m1; (H*-b4]/  
    ex_low_m1=Ex(2); /64jO?mp  
$ tf;\R  
    Ex(KE)=ex_high_m2; VM{`CJ2  
    ex_high_m2=ex_high_m1; "=4`RM  
    ex_high_m1=Ex(KE-1); AH`n  
    %%%%%%%%%%%%%%%%%%计算磁场 ,4y'(DA  
    for i=1:KE-1 9xWC<i  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); F-}-/N]o q  
    end "bZV<;y6  
    plot(Ex); h&4ufx6  
    grid on; PoMkFG6  
    pause(0.01); pT]M]/y/:  
end

我也发到         求助【48小时速问速答区】 » 色散性的FDTD仿真求助  这个里面了，重复发，不算违规吧。

原来是赵博士啊 刚访问了你的科学网博客啊 以后要多多请教

欢迎到这里发帖，好东西是要奖励的

3楼是？ ,Q=)$ `%  
呵呵，不过感兴趣的人好像不多啊。

可不可以说一下Dx和Sx是什么来的呢？ 0q[p{_t`  
最好有相关的物理公式，会更好理解点 #QTfT&m+G}  
谢谢了

要说明的话，可能得整理一份资料出来才行，改天吧。 D:=t*2-Iv  
现在天天写code，发现原来写的真是一个烂字啊。 99<4t$KH  
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); kQ@gO[hS  
这里面竟然不把exp(-1*u0*dt))那些先提取出来，浪费运行时间啊。 9@:BK;Fi  
现在看到乘法反应都是有没有优化的方法，呵呵。

可以先S1与Ex的关系吗？ "u_i[[y  
主要是 C;9t">prk  
     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,%_deV\(  
         Ex(i)=Dx(i)/epsz-Sx1(i); UH[<&v  
         Ex(i)=Ex(i); ^EF

'TO$  
和 j.DHqHx  
    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); sI'a1$  
         Ex(i)=Dx(i)/epsz-Sx1(i); R~)ybf{  
.;9jdGBf  
这两组方程的来历不是太明白，一般看到的等离子波好像都不涉及碰撞速度的，但我在执行的时候发现碰撞速度对结果的影响和等离子振荡频率是差不多的

还有就是步长的设置为什么是dt=dx/2c而不是dt=dx/c呢。但dt=dx/c在这里好像是不稳定的