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

[post]一维等离子体FDTD的Matlab源代码 n,CD4Nv  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s98: *o3  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% D<+ bzC  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% F2Nb5WT  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% v1`bDS?*Q  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% q$e T!'x  
%%%%%%%%%初始化 Z\ "Kd  
clear; Ws2prh^e(  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TKj/6Jz|  
%%%%%%%%%%%%%%%%系统参数 f]@[4<Ny  
u!=]zW%  
TimeT=200;%迭代次数 eyI-s9#t  
KE=2000;%网格树木 m]+X}|  
kc=450;%源的位置 ;`X`c  
kpstart=500;%等离子体开始位置 <) >gg!  
kpstop=1000;%等离子体终止位置 8cYuzt]..  
DispE=zeros(1,TimeT); 34&u]4=L)  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5^G7pI7  
%%%%%%%%%%%%%物理参数 Mqrt-VPh  
c0=3e8;%真空中波速 $ioaunQKP  
f=150e6; C=xo&I7  
f1=75e6; fS;m+D!j@  
lamda=c0/f; X=S}WKu  
WL=50; (E2lv#[  
OMIGA=pi/WL; ?M*C*/R  
zdelta=lamda/WL;%网格大小 `R4W4h'I  
dt=zdelta/(2*c0);%时间间隔 |q1b8A\  
(hD X4
;4  
v.~Nv@+kR  
u0=1e7;%碰撞频率 .#:@cP~v  
fpe=1e7;%等离子体频率 @LwV
mR |{  
wpe=2*pi*fpe;%等离子体圆频率 hW*^1%1  
epsz=1/(4*pi*9*10^9); % 真空介电常数 hr/xpQW  
mu=1/(c0^2*epsz);%磁常数 LH?gJ8`  
ex_low_m1=0; >tE
,8  
ex_low_m2=0; Z_eqM4{  
ex_high_m1=0; s+OvS9et_  
ex_high_m2=0; `Z;B
^Y0  
:O>Nd\UtO  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% YyX^lL_  
%%%%%%%%%%%%%初始化电磁场 IdN%f]=/  
Ex=zeros(1,KE); ~|$) 1  
Ex_Pre=zeros(1,KE); ]ly)z[is"]  
Hy=zeros(1,KE); A;/Xt  
Hy_Pre=zeros(1,KE); QjW~6Z.tI  
Dx=zeros(1,KE); ^^j|0qshL  
Dx_Pre=zeros(1,KE); 7;$L&X  
Sx1=zeros(1,KE); J.CZR[XF#  
Sx2=zeros(1,KE); m\R@.jkZ
 
Sx3=zeros(1,KE); ThT.iD[  
Sx=zeros(1,KE); (_s!,QUe  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% .{-iq(3  
%%%%%%%%%%%%%%%%%%开始计算 Q@3ld6y  
for T=1:TimeT ,=XS%g}l4  
    %%%保存前一时间的电磁场 P~b%;*m}8  
    Ex_Pre=Ex; +E""8kW- Z  
    Hy_Pre=Hy; g*"J10hyP  
    %%%%中间差分计算Dx /_ hfjCE  
    for i=2:KE <E(-QJ  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));  ^qSf  
    end qB`0^V  
    %%%%%%%%加入源 (>)+;$Dr,\  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); {<Xo,U7y  
       -Y!=Iw 4  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); hH|XtQ.n^  
    if T<WL +@<^i?ale  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 SbQ{ >  
    end FrE/K_L  
    %%%计算电场Ex $D2A
in1  
    for i=1:kpstart-1 zL[U;  
        Ex(i)=Dx(i)/epsz; O7L6Htya  
    end MZhJ,km)  
    for i=kpstop+1:KE $7k04e@]  
        Ex(i)=Dx(i)/epsz; ~k:>Xo[|O  
    end 9Rt(G_'  
    if T>50 sVZ}nq{  
        gw=0; 37<GG)  
    end .,iw2:  
    i=kpstart; Q%b46"  
         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); sB*h`vs0T  
         Ex(i)=Dx(i)/epsz-Sx1(i); swe8  
         Ex(i)=Ex(i); Xf'  
     *< SU_dAh  
    for i=kpstart+1:kpstop u.xA}yVS  
         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); ;V
tpq3  
         Ex(i)=Dx(i)/epsz-Sx1(i); 2Xk1AS  
         %Ex(i)=Ex(i); ~jrU#<'G9  
    end k%bTs+]*  
    Sx3=Sx2; 8)2u@sx%  
    Sx2=Sx1; rnt$BB[g  
    Ex(1)=ex_low_m2; aV0;WH_3  
    ex_low_m2=ex_low_m1;  :L+zUlsf  
    ex_low_m1=Ex(2); 36D,el In  
rqG6Ll`=+
 
    Ex(KE)=ex_high_m2; R}=]UOqH-  
    ex_high_m2=ex_high_m1; b "AHw?5F  
    ex_high_m1=Ex(KE-1); s*A|9uf5  
    %%%%%%%%%%%%%%%%%%计算磁场 jak|LOp  
    for i=1:KE-1 < ..
h^3Vd K,  
E'6z7m.  

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% vG=$
UUh@~  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% N1zrfn-VU  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E8V\J  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nWg)zj:  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @G&xq"Fg7  
%%%%%%%%%初始化 |Szr=[  
clear; o".O#^3H%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% `e:RZ
  
%%%%%%%%%%%%%%%%系统参数 s~'C'B?  
TimeT=3000;%迭代次数 )Syf5I  
KE=2000;%网格树木 {@`Uf;hPAX  
kc=450;%源的位置 JN|#  
kpstart=500;%等离子体开始位置 C)dYAq3,8  
kpstop=1000;%等离子体终止位置 o%s}jBo}  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >Qu^{o  
%%%%%%%%%%%%%物理参数 R-0Ohj  
c0=3e8;%真空中波速 tWy<9TF  
zdelta=1e-9;%网格大小 'cCj@bZ9X  
dt=zdelta/(2*c0);%时间间隔 [WSIC *|;  
f=900e12;%Gause脉冲的载频 S~`AnX3!  
d=3e-15%脉冲底座宽度 Or~6t}f  
t0=2.25/f;%脉冲中心时间 v[=E f  
u0=57e12%碰撞频率 Ny<G2!W  
fpe=2000e12;%等离子体频率 L{ ^@O0S  
wpe=2*pi*fpe;%等离子体圆频率 FO3*[O   
epsz=1/(4*pi*9*10^9); % 真空介电常数 n]g,)m  
mu=1/(c0^2*epsz);%磁常数 i2c<q0u  
ex_low_m1=0; 8?R_O}U  
ex_low_m2=0; E^ti!4{<  
ex_high_m1=0; ^d=@RTyo/  
ex_high_m2=0; t|,Ex7  
a0=2*u0/dt+(2/dt)^2; s:.XF|e{  
a1=-8/(dt)^2; eNs
kuG|1  
a2=-2*u0/dt+(2/dt)^2; Q(Y,p`>  
b0=wpe^2+2*u0/dt+(2/dt)^2; F#R\Ot,hv  
b1=2*wpe^2-8/(dt)^2; VZ!$'??  
b2=wpe^2-2*u0/dt+(2/dt)^2; _ q1|\E%`h  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Z-V%lRQ=b  
%%%%%%%%%%%%%初始化电磁场 0(6`dr_  
Ex=zeros(1,KE); ?Gu>!7  
Ex_Pre=zeros(1,KE); xJ.!Q)[  
Hy=zeros(1,KE); y6yseR!  
Hy_Pre=zeros(1,KE); ~#Mx&mZ  
Dx=zeros(1,KE); r\D8_S_  
Dx_Pre=zeros(1,KE); J<O_N~$$*  
Sx1=zeros(1,KE); PA[Rhoit,  
Sx2=zeros(1,KE); -w0>4JDs  
Sx3=zeros(1,KE); M- A}(r +J  
Sx=zeros(1,KE); O/~^}8TLL  
Dx=Ex; JQ4>S<ttJ  
Dx1=Ex; 73<yrBxp  
Dx2=Ex; <08V-   
Dx3=Ex; f@g  
Ex1=Ex; -L3RzX  
Ex2=Ex; }# ^PbM  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BGjTa.
&  
%%%%%%%%%%%%%%%%%%开始计算 LxDhthZi_  
for T=1:TimeT (pg9cM]NA  
    %%%保存前一时间的电磁场 2$UR"P  
    Ex_Pre=Ex; n[-!Jp[  
    Hy_Pre=Hy;
Bbp9Q,4  
    %%%%中间差分计算Dx ,_-*/- 7;8  
    for i=2:KE 1#uw^{n  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); (ytkq(  
    end (`*wiu+i  
    %%%%%%%%加入源 S?tLIi/  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); BS.6d}G4  
     !d)i6W?  
    %%%计算电场Ex Ji=iq=S7  
    for i=1:kpstart-1 m^+~pC5  
        Ex(i)=Dx(i)/epsz; _L.yt5_  
    end gWPa8q<b  
    for i=kpstop+1:KE J8'zvH&I  
        Ex(i)=Dx(i)/epsz; oa7Hx<Y  
    end p+?WhxG)  
    Dx3=Dx2; |g!# \  
Dx2=Dx1; tYa*%|!v  
Dx1=Dx; e8v=n@0  
    for i=kpstart:kpstop %tLq&tyeY  
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); hgltD8,  
    end c^k. <EA  
    Ex2=Ex1; tsD^8~ t|h  
Ex1=Ex; xx8na8  
    Sx3=Sx2; *
@^0xz{\z  
    Sx2=Sx1;  p-k
qX  
    Ex(1)=ex_low_m2; BD+?Ad?  
    ex_low_m2=ex_low_m1; tT:yvU@a  
    ex_low_m1=Ex(2); mws.)  
e:<> Yq+  
    Ex(KE)=ex_high_m2; J>35q'nN]F  
    ex_high_m2=ex_high_m1; qiN'Tuw9  
    ex_high_m1=Ex(KE-1); z(qz(`eGC&  
    %%%%%%%%%%%%%%%%%%计算磁场 2D"/k'iA  
    for i=1:KE-1 }XU- JAn  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); AJu.  
    end !C?z$5g  
    plot(Ex); IF1}}[Ht  
    grid on; "N_?yA#(j  
    pause(0.01); Po3W+;@  
end

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

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

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

3楼是？ q65]bs4M  
呵呵，不过感兴趣的人好像不多啊。

可不可以说一下Dx和Sx是什么来的呢？ BWQ`8
 
最好有相关的物理公式，会更好理解点 Pc=ei  
谢谢了

要说明的话，可能得整理一份资料出来才行，改天吧。 i KQj[%O  
现在天天写code，发现原来写的真是一个烂字啊。 &"JC
8  
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); yQUrHxm  
这里面竟然不把exp(-1*u0*dt))那些先提取出来，浪费运行时间啊。 8^+|I,  
现在看到乘法反应都是有没有优化的方法，呵呵。

可以先S1与Ex的关系吗？ 4IfkYM  
主要是 gM1:*YK  
     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); tQ;Fgv8Y!  
         Ex(i)=Dx(i)/epsz-Sx1(i); st"@kHQ3  
         Ex(i)=Ex(i); KS~Q[-F1P  
和 @]X!#&2>  
    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}7g>  
         Ex(i)=Dx(i)/epsz-Sx1(i);  !vl1#@  
Q_,!(N  
这两组方程的来历不是太明白，一般看到的等离子波好像都不涉及碰撞速度的，但我在执行的时候发现碰撞速度对结果的影响和等离子振荡频率是差不多的

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