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

[post]一维等离子体FDTD的Matlab源代码 s+- aHn  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FJ"9Hs2  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7k|(5P;  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% AoeW<}MO  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% N@)~j+Pz  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BB$>h-M/%#  
%%%%%%%%%初始化 $(=1A>40  
clear; q:^Cw8  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ZRfa!9vl  
%%%%%%%%%%%%%%%%系统参数 yFsXI0I[p  
jRAL(r|  
TimeT=200;%迭代次数 2A+,. S_!x  
KE=2000;%网格树木 *, K \A  
kc=450;%源的位置 xltu g##  
kpstart=500;%等离子体开始位置 {[*_HAy7  
kpstop=1000;%等离子体终止位置 koZp~W-  
DispE=zeros(1,TimeT); ^i\1c-/  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DA~ELje^j  
%%%%%%%%%%%%%物理参数 I_7EfAqg(  
c0=3e8;%真空中波速 wP"|$HN  
f=150e6; F2'cL@E3  
f1=75e6; 9zd/5|W  
lamda=c0/f; /|AuI qW  
WL=50; >~~\==".  
OMIGA=pi/WL; B$EP'5@b  
zdelta=lamda/WL;%网格大小 1A-EP@# J  
dt=zdelta/(2*c0);%时间间隔 ?UDO%`X  
^^uD33@_  
"~ stZ.  
u0=1e7;%碰撞频率 HC,@tfS  
fpe=1e7;%等离子体频率 aS G2K0  
wpe=2*pi*fpe;%等离子体圆频率 F9r*ZyNlx  
epsz=1/(4*pi*9*10^9); % 真空介电常数 P^W47 SO  
mu=1/(c0^2*epsz);%磁常数 c F]3gM  
ex_low_m1=0; Dj"=k
L0  
ex_low_m2=0; :PkZ(WZ9  
ex_high_m1=0; }4%)m  
ex_high_m2=0; ZRGZ'+hw  
y/eX(l<{  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% kH-b!  
%%%%%%%%%%%%%初始化电磁场 3HR]TQ%r  
Ex=zeros(1,KE); p5VSSvV\K  
Ex_Pre=zeros(1,KE); %Oqe7Cx>+  
Hy=zeros(1,KE); #SNI dc>9\  
Hy_Pre=zeros(1,KE); o1+]6s+j}  
Dx=zeros(1,KE); ^?[<!VBI  
Dx_Pre=zeros(1,KE); f om"8iL1  
Sx1=zeros(1,KE); >]8.xkQq  
Sx2=zeros(1,KE); >irT|VTf  
Sx3=zeros(1,KE); rF3]AW(  
Sx=zeros(1,KE); !5t 3Y  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [d>yo_iB  
%%%%%%%%%%%%%%%%%%开始计算 *Bsmn!_cB{  
for T=1:TimeT :Xh`.*{EX  
    %%%保存前一时间的电磁场 rd4'y~#S  
    Ex_Pre=Ex; ,8nZzVo  
    Hy_Pre=Hy; @rE)xco  
    %%%%中间差分计算Dx \RyA}P5S  
    for i=2:KE wJ*
-K-  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); { BL1j  
    end de{YgN  
    %%%%%%%%加入源 tN> B$sv  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); ER1mA:8>E  
       Q.dy $`\  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); N==_'`O1Q0  
    if T<WL ^ZWFj?`\UV  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 }px]  
    end z>;+'>XXgx  
    %%%计算电场Ex MPy][^s!  
    for i=1:kpstart-1 Go7hDmu  
        Ex(i)=Dx(i)/epsz; B/kn&^z$|~  
    end [j U  
    for i=kpstop+1:KE Al=? j#J6p  
        Ex(i)=Dx(i)/epsz; ~U6YN_W  
    end v/.h%6n?  
    if T>50 ?$4R <  
        gw=0; '?3z6%
 
    end (ul-J4E\O  
    i=kpstart; qpqz. {\  
         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); [Fj+p4*N  
         Ex(i)=Dx(i)/epsz-Sx1(i); 8'WMsp
X  
         Ex(i)=Ex(i); q)xl$*g
 
     ;Jn0e:x`E  
    for i=kpstart+1:kpstop ^|i\d\  
         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); mX.3R+t  
         Ex(i)=Dx(i)/epsz-Sx1(i); E816YS='  
         %Ex(i)=Ex(i); yXo0z_ G  
    end G_N-}J>EP  
    Sx3=Sx2; ]O=S2Q  
    Sx2=Sx1; 2*wO5v  
    Ex(1)=ex_low_m2; pQ^,.[[  
    ex_low_m2=ex_low_m1; wW! r}I#  
    ex_low_m1=Ex(2); &
W<>^C2v  
39aCwhh7v  
    Ex(KE)=ex_high_m2; Q>a7Ps@~  
    ex_high_m2=ex_high_m1; <*Kh=v  
    ex_high_m1=Ex(KE-1); 'BdmFKy1  
    %%%%%%%%%%%%%%%%%%计算磁场 eGe[sv"k  
    for i=1:KE-1 < ..  QXxLe*  
Q] yT  

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% m]FaEQVoE  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% {W5ydHXy  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Pz]bZPHn  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3h9Sz8  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n7IL7?!o  
%%%%%%%%%初始化 ~L%Pz0Gg  
clear; oA4D\rn8"  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -}{%Q?rYj  
%%%%%%%%%%%%%%%%系统参数 mje<d"bW  
TimeT=3000;%迭代次数 qwq/Xcv  
KE=2000;%网格树木 (Wd_G-da  
kc=450;%源的位置 
@+'c+
  
kpstart=500;%等离子体开始位置 op hH9D  
kpstop=1000;%等离子体终止位置 Y{B9`Z  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% P?9nTG  
%%%%%%%%%%%%%物理参数 B4]AFRI  
c0=3e8;%真空中波速 <^lJr82  
zdelta=1e-9;%网格大小 -O|&c9W.O  
dt=zdelta/(2*c0);%时间间隔 qbS'|--wH  
f=900e12;%Gause脉冲的载频 v5(q)h  
d=3e-15%脉冲底座宽度 ;i<$7MR.e  
t0=2.25/f;%脉冲中心时间 g%`i=s&N%  
u0=57e12%碰撞频率 ecr886  
fpe=2000e12;%等离子体频率 +>JdYV<?0  
wpe=2*pi*fpe;%等离子体圆频率 ^"3\iA:  
epsz=1/(4*pi*9*10^9); % 真空介电常数 ;~W8v.EW  
mu=1/(c0^2*epsz);%磁常数 L
nP3z5d(  
ex_low_m1=0; aSc{Ft/O  
ex_low_m2=0; drW~)6Lr@  
ex_high_m1=0; kf<c,3A  
ex_high_m2=0; ^vJ"-{  
a0=2*u0/dt+(2/dt)^2; LmP qLH'(Q  
a1=-8/(dt)^2; gks ==|s.  
a2=-2*u0/dt+(2/dt)^2; ,Pn-ZF  
b0=wpe^2+2*u0/dt+(2/dt)^2;
8RVS)D''  
b1=2*wpe^2-8/(dt)^2; ?2{bKIV_  
b2=wpe^2-2*u0/dt+(2/dt)^2; |#{-.r6Y]  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% {jvOHu  
%%%%%%%%%%%%%初始化电磁场 x&'o
]Y  
Ex=zeros(1,KE); KX!i\NHz  
Ex_Pre=zeros(1,KE); k!5m@'f  
Hy=zeros(1,KE); P Q7A~dw9  
Hy_Pre=zeros(1,KE); j5PL
{6  
Dx=zeros(1,KE); &e5(Djz8t  
Dx_Pre=zeros(1,KE); |+>uA[6#  
Sx1=zeros(1,KE); ))!Bg?t-  
Sx2=zeros(1,KE); J , V  
Sx3=zeros(1,KE); n5|l|#c$N  
Sx=zeros(1,KE); J^ `hbP+2  
Dx=Ex; ?AEd(_a!q  
Dx1=Ex; q#MM
  
Dx2=Ex; d')-7C  
Dx3=Ex;  >1A*MP4  
Ex1=Ex; 7KU~(?|:h  
Ex2=Ex; P''X_1oMC  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !=3Rg-'d1  
%%%%%%%%%%%%%%%%%%开始计算 BPqk"HG]T  
for T=1:TimeT  {@gAv!  
    %%%保存前一时间的电磁场 (:5G#?6,  
    Ex_Pre=Ex; ^(ScgoXva  
    Hy_Pre=Hy; 2R]&v;A  
    %%%%中间差分计算Dx s2FngAM;f  
    for i=2:KE 6R`Oh uN.>  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); 'Gjq/L/x  
    end 'n0 .#E_  
    %%%%%%%%加入源 1"}cdq.  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); Ow3P-UzU3  
     #Z\
O}<  
    %%%计算电场Ex  B$^7h!  
    for i=1:kpstart-1 yPV'pT)  
        Ex(i)=Dx(i)/epsz; $6T3y8  
    end FW8-'~  
    for i=kpstop+1:KE Zk&h:c  
        Ex(i)=Dx(i)/epsz; BYi)j6"  
    end 1j0-9Kg'  
    Dx3=Dx2; Q.+|xwz  
Dx2=Dx1; m#f{]+6U  
Dx1=Dx; `2`fiKm  
    for i=kpstart:kpstop &-%X:~|:X  
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); *XbI#L%>  
    end vfcb:x  
    Ex2=Ex1; Tr& }$kird  
Ex1=Ex; N L'R\R  
    Sx3=Sx2; F[saP0 *  
    Sx2=Sx1; '7*=`q{  
    Ex(1)=ex_low_m2; EVt?C+  
    ex_low_m2=ex_low_m1; |%1?3Mpn  
    ex_low_m1=Ex(2); /R
T%0!  
BByCMY  
    Ex(KE)=ex_high_m2; N8v'70  
    ex_high_m2=ex_high_m1; /qU>5;  
    ex_high_m1=Ex(KE-1);
wvMW|  
    %%%%%%%%%%%%%%%%%%计算磁场 y{d^?(-  
    for i=1:KE-1
RKZ6}q1n  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); FB:<zmwR  
    end Zj*\"Ol  
    plot(Ex); GKiukX$'  
    grid on; *1{S*`|cJy  
    pause(0.01); V-!"%fO.s  
end

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

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

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

3楼是？
2RppP?M!  
呵呵，不过感兴趣的人好像不多啊。

可不可以说一下Dx和Sx是什么来的呢？ D-iUN  
最好有相关的物理公式，会更好理解点 "%kGRHq  
谢谢了

要说明的话，可能得整理一份资料出来才行，改天吧。 <r0.ppgY  
现在天天写code，发现原来写的真是一个烂字啊。 uSH>$;a  
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); aX:#'eDB  
这里面竟然不把exp(-1*u0*dt))那些先提取出来，浪费运行时间啊。 U2VV[e)Z!  
现在看到乘法反应都是有没有优化的方法，呵呵。

可以先S1与Ex的关系吗？ #]5&mKi  
主要是 6(ka"Vu~  
     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); _#/!s]$d#  
         Ex(i)=Dx(i)/epsz-Sx1(i); ipx@pNW;"  
         Ex(i)=Ex(i); cmwzKu%  
和 54`bE$:+  
    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); #)W8.  
         Ex(i)=Dx(i)/epsz-Sx1(i); 2neF<H?^o  
XR{5]lKt_  
这两组方程的来历不是太明白，一般看到的等离子波好像都不涉及碰撞速度的，但我在执行的时候发现碰撞速度对结果的影响和等离子振荡频率是差不多的

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