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

原程序就是那个一维PEC边界，正弦波在介质中传播的。 ~K;QdV=YX  
%*********************************************************************** Q7s@,c!m_  
% 1-D FDTD code with PEC boundaries ||Zup\QB  
%*********************************************************************** ^dQ{vL@9b9  
% REUxXaN>Z  
% Program author: Susan C. Hagness L)sgW(@2  
% Department of Electrical and Computer Engineering /'
/I^ab  
% University of Wisconsin-Madison 5>x_G#W  
% 1415 Engineering Drive <g[z jV9p  
% Madison, WI 53706-1691 YT\@fgBt  
% hagness@engr.wisc.edu {a7~P0$  
% .hl_zc#  
% Copyright 2005 rJ!cma  
% ?E([Nc0T  
% This MATLAB M-file implements a finite-difference time-domain Y?0/f[Ax,y  
% solution of Maxwell's curl equations over a one-dimensional space @Wu-&Lb  
% lattice comprised of uniform grid cells. 3A#Tn7  
% & LE5'.s  
% To illustrate the algorithm, a 1-GHz sinusoidal wave propagating i
od%YjZu  
% in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is j_Szw w-  
% modeled. The simplified finite difference system for nonpermeable -pu5O9 @  
% media (discussed in Section 3.6.6 of the text) is implemented. K;?D^n.  
% O`5,L[i1y  
% The grid resolution (dx = 1.5 cm) is chosen to provide 20  ?%,NOX  
% samples per wavelength. The Courant factor S=c*dt/dx is set to [^5\Ww  
% the stability limit (S=1). In 1-D, this is the "magic time step." "j^i6RS  
% The total number of time steps (nmax=240) corresponds to a physical (eN7s_  
% time of 12 ns. C?<pD+]b_  
% jJ,_-ui  
% The grid is terminated with electric-field components at the far-left '0+*  
% (i=1) and far-right (i=ie) boundaries. The sinusoidal wave is >fPo_@O  
% launched by an electric-field hard-source condition at i=1 (see q-F K=r 5  
% Eq. 5.1 in the text). The simplest radiation boundary condition jj5S+ >4  
% for plane wave propagation is used to update the electric field '/W$9jm  
% at i=ie: P49\A^5S!  
% gYKz,$  
% Ez(ie,n+1) = Ez(ie-1,n) nr!N%Hi  
% C 0w+ j  
% To execute this M-file, type "fdtd1D" at the MATLAB prompt. &k}f"TX2  
%
D|;O9iks#  
% This code has been tested in the following Matlab environments: PVCoXOqh  
% Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) U&5*>fd=  
% Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) 2xI|G 3U  
% Matlab version 7.0.0.19920 R14 (May 6, 2004) b9.M'P\  
% Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) Luq4q95]  
% Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) R?{+&r.X  
% /(N/DMl[  
%*********************************************************************** %$!3Pbui  
^J'_CA  
clear .JhQxXj  
FwCb$yE#M  
%*********************************************************************** %ByPwu:f  
% Fundamental constants (`P\nnb  
%*********************************************************************** `9~ %6N?7#  
]?Ef0?44  
cc=2.99792458e8; %speed of light in free space Jup)m/  
muz=4.0*pi*1.0e-7; %permeability of free space h!EA;2yGKa  
epsz=1.0/(cc*cc*muz); %permittivity of free space %q{q.(M#  
'%NglC[J  
freq=1.0e+9; %frequency of source excitation *r7vDc  
lambda=cc/freq; %wavelength of source excitation E[H  
omega=2.0*pi*freq; ',+yD9 @  
G=zWhqieh  
%*********************************************************************** ;a:H-iC  
% Grid parameters >[}oH2oi  
%*********************************************************************** tdy2ZPVtTV  
@AsJnf$y  
ie=200; %number of Ez samples in grid $<}c[Nm  
ih=ie-1; %number of Hy samples in grid 'Uok<;
 
+#W94s~0V  
dx=lambda/20.0; %space increment of 1-D lattice 8K 3dwoT  
dt=dx/cc; %time step (S=1.0) +{J8,^z#  
omegadt=omega*dt; T{YZ`[  
ud1M-lY\U  
nmax=round(12.0e-9/dt); %total number of time steps zzqJeIS  
yK~=6^M  
x=linspace(dx,ih*dx,ih); qL(Q1O!  
<wH+\  
%*********************************************************************** HgI!q<)  
% Material parameters sibYJKOy  
%*********************************************************************** 4fEDg{T  
hp\&g2_S0W  
eps=1.0; g)D_!iz  
sig=5.0e-3; ?^}30V:E  
zw0w."V  
%*********************************************************************** C*7/iRe  
% Updating coefficients for space region with nonpermeable media qNp1<QO0  
%*********************************************************************** k+3qX'fd  
WjV15\,  
scfact=dt/muz/dx; K2   
<xOv8IQ|  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); wQkM:=t5  
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); jc}G+|`  
TJ|Jv8j<s  
%*********************************************************************** UO/sv2CN  
% Field arrays GD W@/oQr  
%*********************************************************************** @mp`C}x"0&  
#(8|9  
ez(1:ie)=0.0; B/? L$m  
hy(1:ih)=0.0; .g*j]!_]  
 6@S6E(^  
%*********************************************************************** [K!9x
M6  
% BEGIN TIME-STEPPING LOOP 'M90Yia  
%*********************************************************************** Tb/TP3N  
8[^'PIz  
for n=1:nmax d0cL9&~qW  
M6Fo.eeK3  
%*********************************************************************** cr7MvXF-  
% Update electric fields eI #Gx_mg  
%*********************************************************************** /7Q|D sa  
x0*{oP  
ez(1)=scfact*sin(omegadt*n); 5j%G7.S\  
G#M)5'Q]U  
rbc=ez(ih); |{jT+  
ez(2:ih)=ca*ez(2:ih)+cb*(hy(2:ih)-hy(1:ih-1)); FR&`
R  
ez(ie)=rbc; n
y={OhP-  
/zn=AAYb  
%*********************************************************************** r~w.J+W  
% Update magnetic fields $o^Z$VmL  
%*********************************************************************** YO6BzS/~  
k9|5TLXq?  
hy(1:ih)=hy(1:ih)+ez(2:ie)-ez(1:ih); [}RoZB&I  
 >@ t  
%*********************************************************************** 9>""x
t
 
% Visualize fields :!ya&o  
%*********************************************************************** t^t% >9o  
Cs,H#L  
rtime=num2str(round(n*dt/1.0e-9)); ZzT
=m*tQ&  
2iAC
_"n  
subplot(2,1,1),plot(x,ez(1:ih)/scfact,'r'),axis([0 3 -1 1]); d_RgKdR )k  
title(['time = ',rtime,' ns']); ylabel('E_z'); I`>U#x*  
subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]); pl[J!d.c  
xlabel('x (meters)');ylabel('H_y'); zM0NRERi  
H=r-f@EOrI  
pause(0.05) >ZA=
9v  
yU$MB,1  
%*********************************************************************** x-^6U  
% END TIME-STEPPING LOOP 4ls:BO;k]  
%*********************************************************************** 2hE(h  
enoj4g7em^  
end WzMYRKZ  
ITu19WG  
我想问的是， p7?  
1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ 1v o)]ff  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 Fe+ @;  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。 iyskAD
S  
有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把 .. 1+o]+Jz|  
"L~(%Nx3  

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

1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ lFBdiIw  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 +X?ErQm  
这两句是吸收电磁波的，所以加了以后应该没有反射波了。 igO>)XbsM  
=u~nLL  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把那个去了，才看到的是没有反射，继续传播啊？ "gJ?LojB<  
对的，去掉以后，边界就相当于电场为0，等同于金属边界了。

另外，我只想简单看真空中的效果等。在原来基础上改动却出现错误。还希望高手帮我看看啊？ (j}7|*.  
或者说，应该怎么改呢？

想看什么样的结果？

我试着回答一下你的这个问题。我没有仔细看，但你可以按我下面说的试试看。 ROFZ*@CH<  
首先，你贴出来的程序是一个很特殊的1D Case，特殊之处在于：空间步长和时间步长取值刚好使得在一个时间步长之内，波形向前传播的距离是一个空间网格。因此程序中省略了很多处理。这是一个带吸收边界条件的FDTD程序。一下对你的问题作答： MtkU]XKGT  
X['9;1Xr  
我想问的是， 0&s6PS%  
1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ '=0}2sF>  
>>这两句是吸收边界条件，不是PEC条件。你后面给出的也是同样一个吸收边界条件。 ;<Q%d~$xy}  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 ~N!HxQ  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。 F
! !HwI
 
有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把那个去了，才看到的是没有反射，继续传播啊？ #gbB
// <  
>>把吸收边界条件去掉，就是PEC条件，也就是强迫最右端的Ez值始终为0。最简单的做法是对最右端的Ez节点不做更新处理。 }JFTe g  
4、另一方面，不去看PEC边界，这几个字。只看编程，我又理解为，把Ez(ie,n+1) = Ez(ie-1,n)这个值赋给边界那个点，所以视觉上应该是？。。。传播过去了。 G
 _cJI  
然后把边界去了，那么Ez(ie)这个点就没有值，就是0，相当于一块金属板，所以反射？ MZG
hN brd  
5、我在这个上，只是想简单看下在真空中的情况，我是把EPS,SIGMA,SCFACT去掉，并且ca=1,cb=dt/epsz,但是出来也不太对。出现错误。这个程序传在附件里了。希望高手给我看看，指出我的问题所在，不甚感激。 Y,s EM%  
>>>>我理解你是将吸收条件和PEC条件搞反了，建议自己编程实现一下，总共也就几十行代码，不要参考你贴出来的那个源码，里面的处理太特殊，容易引起混乱

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

恩，听你的解释，好像我也搞混了。 oI-Fr0!  
是看的上面的PEC边界才接着这么理解的。。。 原来下面的编程和PEC没关系？而是吸收边界的处理。 -}m  
我再回去看看。谢谢了。

我是想看真空中的效果。 P
9\y~W  
我把把EPS,SIGMA,去掉，并且ca=1,cb=scfact*(dt/epsz),但是出来的图形完全不对。 PvUY Q>Kw  
请问问题出在哪呢？ 是作图那边要改么？

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

那个不是相对介质参数的么。 ggTjd"|)  
或者改为 eps=1.0; sigma=0;？貌似这样出来的也不对、 w;$+7  
还有原来的ca，cb不是也需要修改？ qU n>  
ZO1J";>u  
那请教怎么改，是真空下的呢？ U.F65KaKF  
我就是想看到，没有衰减的波，遇到边界反射的情况。谢谢了。