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

原程序就是那个一维PEC边界，正弦波在介质中传播的。 ]:Ocu--  
%*********************************************************************** {#0B~Zr  
% 1-D FDTD code with PEC boundaries sw8Ic
\vT  
%*********************************************************************** enK4`+.7  
% u*}6)=+:  
% Program author: Susan C. Hagness ,5H$Tm,6\S  
% Department of Electrical and Computer Engineering =}u;>[3  
% University of Wisconsin-Madison 4- N>#  
% 1415 Engineering Drive I)O%D3wfMW  
% Madison, WI 53706-1691 c;BQ$je}  
% hagness@engr.wisc.edu Nr6YQH*[  
%
ke'p8Gz  
% Copyright 2005 k;l^wM  
% 9{?<.%  
% This MATLAB M-file implements a finite-difference time-domain ,|?B5n&  
% solution of Maxwell's curl equations over a one-dimensional space OA&NWAm4  
% lattice comprised of uniform grid cells. ^4c2}>f  
% Y-7x**I  
% To illustrate the algorithm, a 1-GHz sinusoidal wave propagating hFs0qPVY  
% in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is DV]Kd 7  
% modeled. The simplified finite difference system for nonpermeable a(BWV?A  
% media (discussed in Section 3.6.6 of the text) is implemented. !V7VM_}@Y  
% ;x|4Tm  
% The grid resolution (dx = 1.5 cm) is chosen to provide 20 GVCyVt[!-  
% samples per wavelength. The Courant factor S=c*dt/dx is set to =KMd! $J\  
% the stability limit (S=1). In 1-D, this is the "magic time step." 
/Y|9!{.  
% The total number of time steps (nmax=240) corresponds to a physical GcHWalm  
% time of 12 ns. l{_>?]S5  
% g}L2\i688  
% The grid is terminated with electric-field components at the far-left j#<#o:If  
% (i=1) and far-right (i=ie) boundaries. The sinusoidal wave is Z:)\j.  
% launched by an electric-field hard-source condition at i=1 (see _f{'&YhUU  
% Eq. 5.1 in the text). The simplest radiation boundary condition GDZe6*  
% for plane wave propagation is used to update the electric field E':y3T@."  
% at i=ie: 4,wdIdSm4  
% 6aXsRhQ~  
% Ez(ie,n+1) = Ez(ie-1,n) H /*^$>0Uo  
% <),FI <~  
% To execute this M-file, type "fdtd1D" at the MATLAB prompt. Q3_ia5 `O  
% ~|R"GloUw  
% This code has been tested in the following Matlab environments: ~$YFfv>  
% Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) `"Tx%>E(U  
% Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) 'R,1Jmx
  
% Matlab version 7.0.0.19920 R14 (May 6, 2004) /,c9&it(M  
% Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) dD@T}^j *|  
% Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) (:vY:-\ bO  
% #oJ9BgDry  
%*********************************************************************** (*^_wq-;  
\Vr(P>  
clear '?wv::t  
1!0BE8s"@  
%*********************************************************************** bmzs!fg_~R  
% Fundamental constants >=r0
94<  
%*********************************************************************** a:H}c9$%  
}eULcgRG  
cc=2.99792458e8; %speed of light in free space 4}LGE>  
muz=4.0*pi*1.0e-7; %permeability of free space ATPc~f  
epsz=1.0/(cc*cc*muz); %permittivity of free space *`=V"nXw$|  
>eWHPO  
freq=1.0e+9; %frequency of source excitation NrhU70y  
lambda=cc/freq; %wavelength of source excitation Gk'J'9*  
omega=2.0*pi*freq; hdfNXZ{A"  
b?h"a<7  
%*********************************************************************** \7jcZ~FBX%  
% Grid parameters 4$1sBY/  
%*********************************************************************** )7*Apy==x  
52{jq18&  
ie=200; %number of Ez samples in grid #eR*|W7o  
ih=ie-1; %number of Hy samples in grid s <Ag8U8  
c R[DT04  
dx=lambda/20.0; %space increment of 1-D lattice AIN_.=]"?  
dt=dx/cc; %time step (S=1.0) kplyZ  
omegadt=omega*dt; b"/P  
7?Xfge%\  
nmax=round(12.0e-9/dt); %total number of time steps y[.0L!C {  
6 @'v6 1'  
x=linspace(dx,ih*dx,ih); & i)p^AmM  

Z\4l+.R`  
%*********************************************************************** s{Ryh.IyI  
% Material parameters ;UArDwH  
%*********************************************************************** i[T!{<  
U3ED3) D  
eps=1.0; OI::0KOv  
sig=5.0e-3; !H~G_?Mf\O  
QU8?/  
%*********************************************************************** Qm^N}>e  
% Updating coefficients for space region with nonpermeable media ){5$8  
%*********************************************************************** n)^B0DnIk  

zG+oZ  
scfact=dt/muz/dx; W29@`93  
oypX.nye_  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); vb\UP&Ip  
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); ,^`+mP  
N=)N   
%*********************************************************************** p@YU7_sF^!  
% Field arrays jbMzcn~ehI  
%*********************************************************************** 7z5AI!s_  
7:9
WiN5b  
ez(1:ie)=0.0; L&.9.Ll  
hy(1:ih)=0.0; 3' mQ=tKa  
WSWaq\9]8  
%*********************************************************************** G0xk @SE  
% BEGIN TIME-STEPPING LOOP *=UEx0_!q  
%*********************************************************************** 7K%Ac  
 G\ru%  
for n=1:nmax lJ:B9n3OzT  
mst-:F[h  
%*********************************************************************** JY4 +MApN  
% Update electric fields ko=v
K%E[  
%***********************************************************************
S^a")U4  
]M-j_("&  
ez(1)=scfact*sin(omegadt*n); x_s9DkX  
Kw"7M~  
rbc=ez(ih); GbBcC#0  
ez(2:ih)=ca*ez(2:ih)+cb*(hy(2:ih)-hy(1:ih-1)); bTb|@  
ez(ie)=rbc; P)7SK&]r;=  
3{]csZvW  
%*********************************************************************** G# .z((Rj  
% Update magnetic fields g.i
iT/b  
%*********************************************************************** k`'^e
/  
SHIK=&\~-  
hy(1:ih)=hy(1:ih)+ez(2:ie)-ez(1:ih); u01x}Ff~6  
W7_X=>l  
%*********************************************************************** j^Bo0{{  
% Visualize fields WVir[Kv%  
%*********************************************************************** yAW%y  
_ ?xORzO  
rtime=num2str(round(n*dt/1.0e-9)); <t.yn\G-w  
h/QZcA  
subplot(2,1,1),plot(x,ez(1:ih)/scfact,'r'),axis([0 3 -1 1]); EO:i+e]=  
title(['time = ',rtime,' ns']); ylabel('E_z'); C8e{9CF  
subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]); !J3g,p*  
xlabel('x (meters)');ylabel('H_y'); diaLw  
n]r7} 2hM  
pause(0.05) }BzV<8F  
FI Io{ru  
%*********************************************************************** 7??+8T#n*  
% END TIME-STEPPING LOOP :-69,e  
%*********************************************************************** l]uF!']f  
5)AMl)  
end }c:s+P+/  
jLM1~`&  
我想问的是， Vyf r>pgW1  
1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ VNp[J'a>VZ  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 ;S
agN  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。 K\,)9:`t  
有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把 .. Nw/4z$].J  
'[I?G6  

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

1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ Kh&W\\K  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 .);~H#  
这两句是吸收电磁波的，所以加了以后应该没有反射波了。 ndg1E;>  
gR1vUad7  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把那个去了，才看到的是没有反射，继续传播啊？  bi/ AQ^  
对的，去掉以后，边界就相当于电场为0，等同于金属边界了。

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

想看什么样的结果？

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

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

恩，听你的解释，好像我也搞混了。 &> _aY
#  
是看的上面的PEC边界才接着这么理解的。。。 原来下面的编程和PEC没关系？而是吸收边界的处理。 FRL;fF  
我再回去看看。谢谢了。

我是想看真空中的效果。 RQu[FZT,  
我把把EPS,SIGMA,去掉，并且ca=1,cb=scfact*(dt/epsz),但是出来的图形完全不对。 8
M,z#DF  
请问问题出在哪呢？ 是作图那边要改么？

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

那个不是相对介质参数的么。 gO@LJ  
或者改为 eps=1.0; sigma=0;？貌似这样出来的也不对、 M6V^ur 1  
还有原来的ca，cb不是也需要修改？ dYlVJ_0Zr  
?+%bEZ`  
那请教怎么改，是真空下的呢？ q$`>[&I~)  
我就是想看到，没有衰减的波，遇到边界反射的情况。谢谢了。