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

[post]一维等离子体FDTD的Matlab源代码 1$cl "d`~  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% VuLb9Kn  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Qt u;
_  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T CT8OU|  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rfV'EjiM}  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5Fydh0.  
%%%%%%%%%初始化  %:26v  
clear; M+"6VtZH  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% <Bo\a3Z  
%%%%%%%%%%%%%%%%系统参数 a)|y0w)vV  
N^ +q^iW  
TimeT=200;%迭代次数 2gWR2 H@  
KE=2000;%网格树木 Mo/R+\u+Y  
kc=450;%源的位置 DJGafX^  
kpstart=500;%等离子体开始位置 nDi^s{  
kpstop=1000;%等离子体终止位置 !ooi.Oz*Tu  
DispE=zeros(1,TimeT); 2EgvS!"  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ~EtGR # N  
%%%%%%%%%%%%%物理参数 kBP?_ O  
c0=3e8;%真空中波速 ?K$&|w%{3  
f=150e6; s|IBX0^@  
f1=75e6; Y9BQLu4F  
lamda=c0/f; Om.%K>V  
WL=50; &5 7c!)  
OMIGA=pi/WL; B*/!s7c.  
zdelta=lamda/WL;%网格大小 9J:|"@)N  
dt=zdelta/(2*c0);%时间间隔 b'wy{~l@  
7_~sa{1R.  
x,nl PU  
u0=1e7;%碰撞频率  \? /'
 
fpe=1e7;%等离子体频率 <B@NSj  
wpe=2*pi*fpe;%等离子体圆频率 $(}rTm  
epsz=1/(4*pi*9*10^9); % 真空介电常数 .2I?^w&j+  
mu=1/(c0^2*epsz);%磁常数 r8"2C#  
ex_low_m1=0; S1|5+PPs  
ex_low_m2=0; ]p|?S[!=  
ex_high_m1=0; |s3;`Nxu7  
ex_high_m2=0; .!q_jl%U  
u|KjoO   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Na@bXcz)  
%%%%%%%%%%%%%初始化电磁场 'Agw~ &$  
Ex=zeros(1,KE); jCY
~Wc  
Ex_Pre=zeros(1,KE); ,
#;hI{E  
Hy=zeros(1,KE); !mv5i%3  
Hy_Pre=zeros(1,KE); Nuqmp7C  
Dx=zeros(1,KE); tE%g)hL-  
Dx_Pre=zeros(1,KE); e5mu-  
Sx1=zeros(1,KE); 1zRYd`IPoq  
Sx2=zeros(1,KE); 2n.HmS  
Sx3=zeros(1,KE); ~\=D@G,9  
Sx=zeros(1,KE); Zk`y"[J  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NV5qF/<M  
%%%%%%%%%%%%%%%%%%开始计算 5V 2ZAYV  
for T=1:TimeT A%#M#hD/  
    %%%保存前一时间的电磁场 I`{3I-E  
    Ex_Pre=Ex; -!!]1\S*Y  
    Hy_Pre=Hy; 9!AvsC9  
    %%%%中间差分计算Dx
S(@kdL  
    for i=2:KE U\%r33L )  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); 2l?^\9&
 
    end G=y~)B}  
    %%%%%%%%加入源 SM~~:  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); ")KqPD6k  
       h6k" D4o\  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); 3$hIc)  
    if T<WL )\yK61aX  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 YCRE-5!  
    end A=kOSq 4Q  
    %%%计算电场Ex `E|i
8M3g  
    for i=1:kpstart-1 2hV -h  
        Ex(i)=Dx(i)/epsz; 'p5M|h\:T  
    end J0V m&TY  
    for i=kpstop+1:KE a&{Y~Og?%  
        Ex(i)=Dx(i)/epsz; 6uD<E  
    end 1 b7jNkQ  
    if T>50 :]:)c8!6  
        gw=0; iw#~xel<ez  
    end JuZkE9C,${  
    i=kpstart; ;PaU"z+Je~  
         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); H!*ypJ  
         Ex(i)=Dx(i)/epsz-Sx1(i); 0SvPr[ >  
         Ex(i)=Ex(i); s[GHDQ;!  
     Oj-\  
    for i=kpstart+1:kpstop +iQ@J+k  
         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); .I_atv  
         Ex(i)=Dx(i)/epsz-Sx1(i); :G>w MMv&z  
         %Ex(i)=Ex(i); brp3xgQ`]  
    end s(py7{ ^K  
    Sx3=Sx2; YM`T"`f  
    Sx2=Sx1; E5i5gE"\  
    Ex(1)=ex_low_m2; <ll?rPio"  
    ex_low_m2=ex_low_m1; P;"moluE;  
    ex_low_m1=Ex(2); mr7Oi `dE  
+TbAtkEF*  
    Ex(KE)=ex_high_m2; xHt7/8wF  
    ex_high_m2=ex_high_m1; 1j<uFhi>  
    ex_high_m1=Ex(KE-1); YE@yts  
    %%%%%%%%%%%%%%%%%%计算磁场 NsI.mTc2  
    for i=1:KE-1 < .. 0{vT`e'  
uZ8-?  

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ~588M 8~  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2Snb+,o2  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% .KKecdd?=  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ^7\kvW  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Uk]jy>7;!  
%%%%%%%%%初始化 ^85Eveu  
clear; ~\(c;J*Ir  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Eo2`Vr9g  
%%%%%%%%%%%%%%%%系统参数 R_B0CM<!  
TimeT=3000;%迭代次数 @=l6zd@  
KE=2000;%网格树木 l,lqhq\  
kc=450;%源的位置 _|US`,kfc  
kpstart=500;%等离子体开始位置 %JrZ
Ms>  
kpstop=1000;%等离子体终止位置 6&0@k^7~  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s8&q8r7%  
%%%%%%%%%%%%%物理参数 !O%!A<3  
c0=3e8;%真空中波速 R/x3+_.f  
zdelta=1e-9;%网格大小 d7](fw@c  
dt=zdelta/(2*c0);%时间间隔 d)1gpRp  
f=900e12;%Gause脉冲的载频 !*Is0``  
d=3e-15%脉冲底座宽度 Vm<_
e  
t0=2.25/f;%脉冲中心时间 F/ZFO5C%  
u0=57e12%碰撞频率 <&Xl b0  
fpe=2000e12;%等离子体频率 B>c$AS\5y  
wpe=2*pi*fpe;%等离子体圆频率 l,hOnpm9  
epsz=1/(4*pi*9*10^9); % 真空介电常数 UH-873AK  
mu=1/(c0^2*epsz);%磁常数 l#enbQ`-~  
ex_low_m1=0; -9FGFBm4]  
ex_low_m2=0; `$Rgn3  
ex_high_m1=0; *VhEl7  
ex_high_m2=0; 5c3-?u!  
a0=2*u0/dt+(2/dt)^2; ?'0!>EjY"  
a1=-8/(dt)^2; X PyDZk/m  
a2=-2*u0/dt+(2/dt)^2; LUD.  
b0=wpe^2+2*u0/dt+(2/dt)^2; !DOyOTR&3  
b1=2*wpe^2-8/(dt)^2; z |llf7:  
b2=wpe^2-2*u0/dt+(2/dt)^2; 8)bR\s  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% pWKE`x^  
%%%%%%%%%%%%%初始化电磁场 2c]"*Pb  
Ex=zeros(1,KE); Q&.uL}R  
Ex_Pre=zeros(1,KE); "7y,d%H  
Hy=zeros(1,KE); 2|^@=.4\  
Hy_Pre=zeros(1,KE); T='uqKW\  
Dx=zeros(1,KE); H;4QuB'^  
Dx_Pre=zeros(1,KE); mq[=,,#  
Sx1=zeros(1,KE); iGSJ\  
Sx2=zeros(1,KE); OM0r*<D"!  
Sx3=zeros(1,KE); uCr& `  
Sx=zeros(1,KE); K7`6G[RMb  
Dx=Ex; `sqr>QD  
Dx1=Ex; %Zk6K!MY#  
Dx2=Ex; j yD3Sa3  
Dx3=Ex; LH2B*8=^2  
Ex1=Ex; 1l$C3c  
Ex2=Ex; D%]S>g5
k  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% AwNr}9`  
%%%%%%%%%%%%%%%%%%开始计算 ]uox ^HC  
for T=1:TimeT Qpv#&nfUi6  
    %%%保存前一时间的电磁场 z(LR!hr  
    Ex_Pre=Ex; '27$x&6>S  
    Hy_Pre=Hy; GIzB1cl:  
    %%%%中间差分计算Dx uQ-GJI^t  
    for i=2:KE vQLYWRXiA  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); <z
\SKR[  
    end YA$YT8iMe  
    %%%%%%%%加入源 6}-No  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); R?iCJ5m  
     86)2\uan  
    %%%计算电场Ex *"WP*A\1  
    for i=1:kpstart-1 @ &
N  
        Ex(i)=Dx(i)/epsz; P4Pc;8T@!  
    end Q~nVbj?c2v  
    for i=kpstop+1:KE g6%]uCFB  
        Ex(i)=Dx(i)/epsz; IMwV9rF  
    end nQmHYOF%  
    Dx3=Dx2; 'Wnh1|z  
Dx2=Dx1; RJ@79L*#  
Dx1=Dx; nRc\!4  
    for i=kpstart:kpstop Cd)g8
<  
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); } doAeTZ  
    end I :<,9.   
    Ex2=Ex1; pFS@yHs  
Ex1=Ex; sMG
o1pG(  
    Sx3=Sx2; -$<oY88  
    Sx2=Sx1; _aevaWtEx  
    Ex(1)=ex_low_m2; Y M:9m)  
    ex_low_m2=ex_low_m1; R
b:H3zh  
    ex_low_m1=Ex(2); BS fmS(.  
5NZuaN  
    Ex(KE)=ex_high_m2; Jm<NDE~rw  
    ex_high_m2=ex_high_m1; IOZw[9](+  
    ex_high_m1=Ex(KE-1); ^g*Sy, A  
    %%%%%%%%%%%%%%%%%%计算磁场 C33Jzn's  
    for i=1:KE-1 LH(P<k&  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); ybiTWM  
    end AB/${RGf+  
    plot(Ex); aC[G_ACwc  
    grid on; J[:#(c&c!1  
    pause(0.01); k)-+ZmMOh  
end

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

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

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

3楼是？ <[mvfw  
呵呵，不过感兴趣的人好像不多啊。

可不可以说一下Dx和Sx是什么来的呢？ u-a*fT  
最好有相关的物理公式，会更好理解点 mGmkeD'  
谢谢了

要说明的话，可能得整理一份资料出来才行，改天吧。 Mb6#97  
现在天天写code，发现原来写的真是一个烂字啊。 D2`tWRm0  
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?a2P6mAh  
这里面竟然不把exp(-1*u0*dt))那些先提取出来，浪费运行时间啊。 ydCVG,"  
现在看到乘法反应都是有没有优化的方法，呵呵。

可以先S1与Ex的关系吗？ ~c~$2Xo  
主要是 _pSCv:3T  
     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); M{U{iS  
         Ex(i)=Dx(i)/epsz-Sx1(i); cTO\
Vhg  
         Ex(i)=Ex(i); d,fX3  
和 sh []OSM
 
    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); ,6,sz]3-  
         Ex(i)=Dx(i)/epsz-Sx1(i); O2|[g8(_F  
-c-#1_X5  
这两组方程的来历不是太明白，一般看到的等离子波好像都不涉及碰撞速度的，但我在执行的时候发现碰撞速度对结果的影响和等离子振荡频率是差不多的

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