[求助]一维Matlab仿真的几个问题

原程序就是那个一维PEC边界，正弦波在介质中传播的。 Wk7E&?-:6  
%*********************************************************************** L3HC-  
% 1-D FDTD code with PEC boundaries t O.5  
%*********************************************************************** !AJkd.  
% -5  
% Program author: Susan C. Hagness Bk3\NPa  
% Department of Electrical and Computer Engineering n=4  
% University of Wisconsin-Madison [/q Bvuun  
% 1415 Engineering Drive 0
ZwXuq  
% Madison, WI 53706-1691 T0dD:sN  
% hagness@engr.wisc.edu bwhH2^ !  
% 2DPv7\fW  
% Copyright 2005 jZ-s6r2=  
% q/zU'7%@  
% This MATLAB M-file implements a finite-difference time-domain m W>Iib|  
% solution of Maxwell's curl equations over a one-dimensional space aL[6}U0(}  
% lattice comprised of uniform grid cells. 6\I^]\YO  
% w!H(zjv&(  
% To illustrate the algorithm, a 1-GHz sinusoidal wave propagating ?2#'>B  
% in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is czIAx1R9  
% modeled. The simplified finite difference system for nonpermeable wdP(MkaV  
% media (discussed in Section 3.6.6 of the text) is implemented. \~A qA!)6  
% 6d/Q"As  
% The grid resolution (dx = 1.5 cm) is chosen to provide 20 wH!$TAZ:Yw  
% samples per wavelength. The Courant factor S=c*dt/dx is set to b MD|  
% the stability limit (S=1). In 1-D, this is the "magic time step." "G%</G8M  
% The total number of time steps (nmax=240) corresponds to a physical ssLswb  
% time of 12 ns. 2#:p:R8I>  
% 5b/ ~]v  
% The grid is terminated with electric-field components at the far-left .B<Bqr@?8  
% (i=1) and far-right (i=ie) boundaries. The sinusoidal wave is E\DA3lq  
% launched by an electric-field hard-source condition at i=1 (see 7^#f)Vp  
% Eq. 5.1 in the text). The simplest radiation boundary condition 1xEOYM)  
% for plane wave propagation is used to update the electric field 4@{?4k-cq  
% at i=ie: X-nC2[tu'W  
% 9MfU{
4:;I  
% Ez(ie,n+1) = Ez(ie-1,n) OWwqCPz.  
% Z6${nUX  
% To execute this M-file, type "fdtd1D" at the MATLAB prompt. d2Q*1Q@u  
% 8cOft ;|qB  
% This code has been tested in the following Matlab environments: @I-gs(  
% Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) %H\J@{f  
% Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) [5~mP`He  
% Matlab version 7.0.0.19920 R14 (May 6, 2004) -_Z4)"k  
% Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) Wgh@XB  
% Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) ;$&\:-6A#  
% S=-$:65  
%*********************************************************************** p&RC#wYu  
;<Z6Y3>I8  
clear H}kSXKO8!8  
>nSt<e  
%*********************************************************************** Kw`CN  
% Fundamental constants qJ`:$U  
%*********************************************************************** `K
5*Fjx  
{*B0lr`  
cc=2.99792458e8; %speed of light in free space xrvM}Il  
muz=4.0*pi*1.0e-7; %permeability of free space 4zS0kk;+  
epsz=1.0/(cc*cc*muz); %permittivity of free space Y&yfm/Ru  
\AroSy9  
freq=1.0e+9; %frequency of source excitation ,\v'%,:C  
lambda=cc/freq; %wavelength of source excitation 3E*m.jX  
omega=2.0*pi*freq; Jf?6y~X>Y  
2lsUCQI;  
%*********************************************************************** R6(:l; W  
% Grid parameters Gqd|F>  
%*********************************************************************** Bz_'>6w  
~5&4
s  
ie=200; %number of Ez samples in grid i:aW .QZ.  
ih=ie-1; %number of Hy samples in grid uEJ8Lmi
 
:sg}e  
dx=lambda/20.0; %space increment of 1-D lattice seEo)m`d  
dt=dx/cc; %time step (S=1.0) T%)E!:}v  
omegadt=omega*dt; k2v:F  
_

xv3UzD  
nmax=round(12.0e-9/dt); %total number of time steps ?< b{  
ecs 0iW-,  
x=linspace(dx,ih*dx,ih); !\4B.  
_Z[0:4  
%*********************************************************************** GqRXNs!  
% Material parameters <EUR:  
%*********************************************************************** 9r]|P}yuS  
mF\!~ag|  
eps=1.0; ?MRY*[$  
sig=5.0e-3; #{|cSaX<  
70 7( LG  
%*********************************************************************** ;9OhK71}  
% Updating coefficients for space region with nonpermeable media Tp;W4]'a*:  
%*********************************************************************** cw!,.o%cD  
Oh$:qu7o0&  
scfact=dt/muz/dx; WuU
wd#e  
c$ZVvu  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); =^u;uS[IW  
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); ??p%_{QY~b  
@p*)^D6E\  
%*********************************************************************** To>,8E+GAb  
% Field arrays GAgTy  
%***********************************************************************
a,vS{434J  
0uDDaFS  
ez(1:ie)=0.0; lcZ.}   
hy(1:ih)=0.0; 2.z-&lFBZ  
;WSW&2  
%*********************************************************************** zw+aZDcV(  
% BEGIN TIME-STEPPING LOOP >E+g.5 ,:W  
%*********************************************************************** L_!ShE
  
h[]9F.[  
for n=1:nmax RJ3oI+gI  
0mSP  
%*********************************************************************** '3672wF/  
% Update electric fields ;Z{D@g+  
%*********************************************************************** 4`#Q  
QnVr)4"  
ez(1)=scfact*sin(omegadt*n); 2g{tzR_j  
\_1a#|97e  
rbc=ez(ih); x>[]Qk^?q  
ez(2:ih)=ca*ez(2:ih)+cb*(hy(2:ih)-hy(1:ih-1)); %BGg?&  
ez(ie)=rbc; !}TsFa  
y,nmPX?]n  
%*********************************************************************** |2q3spd  
% Update magnetic fields AVpg  
%*********************************************************************** \Vf:/9^  
2HFn\kjj.s  
hy(1:ih)=hy(1:ih)+ez(2:ie)-ez(1:ih); D|9+:Y  
12n:)yQy  
%*********************************************************************** jCJcVO>OZ  
% Visualize fields 4MS<t FH)  
%*********************************************************************** +\Vm t[v  
h`|04Q  
rtime=num2str(round(n*dt/1.0e-9)); 3{3@>8{w  
w95M B*N  
subplot(2,1,1),plot(x,ez(1:ih)/scfact,'r'),axis([0 3 -1 1]); ?+d`_/IB  
title(['time = ',rtime,' ns']); ylabel('E_z'); !CYC7HeF  
subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]); \
2s`mCY  
xlabel('x (meters)');ylabel('H_y'); 3^y(@XFt  
?A3L8^tR  
pause(0.05) !e|\1v'
0  
Eu?z!  
%*********************************************************************** pIlEoG=[_  
% END TIME-STEPPING LOOP )SmnLvL  
%*********************************************************************** p' >i3T(  
Q|AZv>'!  
end &|>~7(
 
5X!-Hj  
我想问的是， ZkbE&7Z  
1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ lv 8EfN  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 SL4?E<J
b  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。 qG6s.TcG  
有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把 .. f,HUr% @  
iT-coI  

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

1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ YFO{i-*q  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 `geHSx_  
这两句是吸收电磁波的，所以加了以后应该没有反射波了。 ":Wq<Z'  
oh9
 ;_~  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把那个去了，才看到的是没有反射，继续传播啊？ #4^d#Gj  
对的，去掉以后，边界就相当于电场为0，等同于金属边界了。

另外，我只想简单看真空中的效果等。在原来基础上改动却出现错误。还希望高手帮我看看啊？ G95,J/w  
或者说，应该怎么改呢？

想看什么样的结果？

我试着回答一下你的这个问题。我没有仔细看，但你可以按我下面说的试试看。 hlTbCl  
首先，你贴出来的程序是一个很特殊的1D Case，特殊之处在于：空间步长和时间步长取值刚好使得在一个时间步长之内，波形向前传播的距离是一个空间网格。因此程序中省略了很多处理。这是一个带吸收边界条件的FDTD程序。一下对你的问题作答： 6_LeP9s )  
bS.w<V Ew  
我想问的是， Wboh2:TH:  
1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ Ucj?$=  
>>这两句是吸收边界条件，不是PEC条件。你后面给出的也是同样一个吸收边界条件。 niVR!l  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 5E:$\z;  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。 :|7#D,2  
有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把那个去了，才看到的是没有反射，继续传播啊？ v9$!v^U"D  
>>把吸收边界条件去掉，就是PEC条件，也就是强迫最右端的Ez值始终为0。最简单的做法是对最右端的Ez节点不做更新处理。 {h<D/:^v  
4、另一方面，不去看PEC边界，这几个字。只看编程，我又理解为，把Ez(ie,n+1) = Ez(ie-1,n)这个值赋给边界那个点，所以视觉上应该是？。。。传播过去了。 <=uYfi3,  
然后把边界去了，那么Ez(ie)这个点就没有值，就是0，相当于一块金属板，所以反射？ K|& f5w  
5、我在这个上，只是想简单看下在真空中的情况，我是把EPS,SIGMA,SCFACT去掉，并且ca=1,cb=dt/epsz,但是出来也不太对。出现错误。这个程序传在附件里了。希望高手给我看看，指出我的问题所在，不甚感激。 v*;d  
>>>>我理解你是将吸收条件和PEC条件搞反了，建议自己编程实现一下，总共也就几十行代码，不要参考你贴出来的那个源码，里面的处理太特殊，容易引起混乱

最简单的，比如模拟一维真空中波的传播。只是看看而已，出不来才奇怪。呵呵。

恩，听你的解释，好像我也搞混了。 dBWny&  
是看的上面的PEC边界才接着这么理解的。。。 原来下面的编程和PEC没关系？而是吸收边界的处理。 Z
9{~t  
我再回去看看。谢谢了。

我是想看真空中的效果。 U%3N=M  
我把把EPS,SIGMA,去掉，并且ca=1,cb=scfact*(dt/epsz),但是出来的图形完全不对。 m^X51,+<  
请问问题出在哪呢？ 是作图那边要改么？

为什么要把eps,sigma去掉呢？

那个不是相对介质参数的么。 l7`{O/hN  
或者改为 eps=1.0; sigma=0;？貌似这样出来的也不对、 Jn+-G4h$  
还有原来的ca，cb不是也需要修改？ x`E<]z*w}  
sx?IIF
F  
那请教怎么改，是真空下的呢？ }rQQe:{]B  
我就是想看到，没有衰减的波，遇到边界反射的情况。谢谢了。