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

原程序就是那个一维PEC边界，正弦波在介质中传播的。 AD4L`0D  
%*********************************************************************** e*'|iuDrY  
% 1-D FDTD code with PEC boundaries 4jyr\=42F'  
%*********************************************************************** JWm^RQ  
% JQVw6*u{  
% Program author: Susan C. Hagness ]oWZ{#r2  
% Department of Electrical and Computer Engineering | 9\7xT  
% University of Wisconsin-Madison cM7k){  
% 1415 Engineering Drive 7! A%6  
% Madison, WI 53706-1691 `__?7"p )\  
% hagness@engr.wisc.edu Eg-Mm4o  
% BD-c 0-+m  
% Copyright 2005 .]sIoB-54  
% B$[%pm`'2  
% This MATLAB M-file implements a finite-difference time-domain R*TGn_J`  
% solution of Maxwell's curl equations over a one-dimensional space R4rm>zisVX  
% lattice comprised of uniform grid cells. {dr&46$p  
% <:yq~?  
% To illustrate the algorithm, a 1-GHz sinusoidal wave propagating  lN`_0  
% in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is 2ZzD^:V[}  
% modeled. The simplified finite difference system for nonpermeable ~)_
?:.Da  
% media (discussed in Section 3.6.6 of the text) is implemented. BC0c c[x  
% 6/WK((Fd  
% The grid resolution (dx = 1.5 cm) is chosen to provide 20 '|A5a+[  
% samples per wavelength. The Courant factor S=c*dt/dx is set to Gk]qE]hi  
% the stability limit (S=1). In 1-D, this is the "magic time step." FsPDWy&x  
% The total number of time steps (nmax=240) corresponds to a physical _)Z7Le:f!  
% time of 12 ns. $rQFM[  
% :"+UG-S$6  
% The grid is terminated with electric-field components at the far-left ]OCJ~Zw  
% (i=1) and far-right (i=ie) boundaries. The sinusoidal wave is J7xT6Q=  
% launched by an electric-field hard-source condition at i=1 (see ("M#R!3  
% Eq. 5.1 in the text). The simplest radiation boundary condition 8UY=}R2C
  
% for plane wave propagation is used to update the electric field n4_:#L?  
% at i=ie: Ev|{~U  
% EwBN+v;)  
% Ez(ie,n+1) = Ez(ie-1,n) tP^mq>  
% p31rhe   
% To execute this M-file, type "fdtd1D" at the MATLAB prompt. g KmRjK  
% R_*D7|v  
% This code has been tested in the following Matlab environments: j?KB8oY`TP  
% Matlab version 6.1.0.450 Release 12.1 (May 18, 2001)
$?JLCa  
% Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) 'V9aB
5O&  
% Matlab version 7.0.0.19920 R14 (May 6, 2004) E<G@LT  
% Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) a]=vq(N'r  
% Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) "i<3}6/*  
% -O>mY)  
%*********************************************************************** "drh+oo.  
`zOAltfd  
clear IWRq:Gw  
n#L2cv~Aj"  
%*********************************************************************** SUi1*S  
% Fundamental constants ~m09yc d<  
%*********************************************************************** /B
?SaKh  
j6d"8oH _  
cc=2.99792458e8; %speed of light in free space byj mH  
muz=4.0*pi*1.0e-7; %permeability of free space V-U  ^O45  
epsz=1.0/(cc*cc*muz); %permittivity of free space /$.vHt5nt  
w wRT$-!  
freq=1.0e+9; %frequency of source excitation SoGLsO+R  
lambda=cc/freq; %wavelength of source excitation Rc.<0#  
omega=2.0*pi*freq; _x|8U'|Ce  
]!
J3?G  
%*********************************************************************** jl0Eg  
% Grid parameters 'GdlqbX(%  
%*********************************************************************** {F9Qy0.*u  
[ *a>{sO[  
ie=200; %number of Ez samples in grid Nb9V/2c;V  
ih=ie-1; %number of Hy samples in grid 4Z p5o`*g2  
]\mb6Hc  
dx=lambda/20.0; %space increment of 1-D lattice wj5s5dH  
dt=dx/cc; %time step (S=1.0) ,4T$  
omegadt=omega*dt; qsRfG~Cg  
yc4f\0B/  
nmax=round(12.0e-9/dt); %total number of time steps kkBV
;v%a  
Ac(irPrD  
x=linspace(dx,ih*dx,ih); LW 3J$Am  
|3lAye,t)a  
%*********************************************************************** D}/.;]w<[&  
% Material parameters >,]e[/p  
%*********************************************************************** -/7=\kao%  
VGUDUM.8  
eps=1.0; 
>sS:x,-  
sig=5.0e-3; ?$&rC0t  
it|:P  
%*********************************************************************** JOne&{h]J"  
% Updating coefficients for space region with nonpermeable media -:p1gg&  
%*********************************************************************** 9X&qdA/q  

f I-"8f0_  
scfact=dt/muz/dx; j#>![km Mu  
sU_4+Mk  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); whZ],R*u  
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); h)RM9813<  
ao5yW;^y  
%*********************************************************************** c1!/jTX$  
% Field arrays 4k?JxA)  
%*********************************************************************** <mgTWv  
N$a-i  
ez(1:ie)=0.0; >Pd23TsN  
hy(1:ih)=0.0; @.@#WHde  
ITqigGan%  
%*********************************************************************** *rKv`nva5  
% BEGIN TIME-STEPPING LOOP Z q>.;>  
%*********************************************************************** eKti+n.  
5@rqU(]<  
for n=1:nmax y\|\9Q%D  
WK)k-A^q  
%*********************************************************************** |nZB/YZt  
% Update electric fields
4qcIoO  
%*********************************************************************** AS"|r  
VRP.tD  
ez(1)=scfact*sin(omegadt*n); @O0vh$3t0  
R/xCS.yl}  
rbc=ez(ih); 'xI+kyu  
ez(2:ih)=ca*ez(2:ih)+cb*(hy(2:ih)-hy(1:ih-1)); Uk ;.Hrt.  
ez(ie)=rbc; [a*>@IR  
Kf<_A{s  
%*********************************************************************** #_sVB~sn@  
% Update magnetic fields 1% %Tm"  
%*********************************************************************** <).qe Z  
4xn^`xf9  
hy(1:ih)=hy(1:ih)+ez(2:ie)-ez(1:ih); /,%o<Ql9  
aNU%OeQA  
%*********************************************************************** HZ4 ^T7G  
% Visualize fields 'wq:F?viF  
%*********************************************************************** ^M5uLm-_s  
R*G>)YH  
rtime=num2str(round(n*dt/1.0e-9)); 0(f;am0y  
ly5L-=Xb  
subplot(2,1,1),plot(x,ez(1:ih)/scfact,'r'),axis([0 3 -1 1]); \~j(ui|  
title(['time = ',rtime,' ns']); ylabel('E_z'); 'eY[?LJ]U  
subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]); PCnJ2  
xlabel('x (meters)');ylabel('H_y'); X1qj l_A  

Y6;9j=[  
pause(0.05) \gv x)S11  
/j\TmcnU^  
%*********************************************************************** &L`p4AZ  
% END TIME-STEPPING LOOP N!&VBx^z  
%*********************************************************************** yVv3S[J  
$K6`Q4`  
end SFNd,(kB*z  
 #X_M  
我想问的是， Xn'>k[}<k  
1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ $^ dk>Hj>4  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 / hdl  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。 k6XmBBIj-  
有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把 .. <Py/uF|  
(\Zo"x;(  

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

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

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

想看什么样的结果？

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

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

恩，听你的解释，好像我也搞混了。 LR^b
?.#>  
是看的上面的PEC边界才接着这么理解的。。。 原来下面的编程和PEC没关系？而是吸收边界的处理。 }vL[N~5\  
我再回去看看。谢谢了。

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

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

那个不是相对介质参数的么。 HHT K{X+  
或者改为 eps=1.0; sigma=0;？貌似这样出来的也不对、 +|\dVe.  
还有原来的ca，cb不是也需要修改？ wS8qua  
:h?Zg(l  
那请教怎么改，是真空下的呢？ C <]rY  
我就是想看到，没有衰减的波，遇到边界反射的情况。谢谢了。