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

原程序就是那个一维PEC边界，正弦波在介质中传播的。 ALieUf  
%*********************************************************************** IqJ=\  
% 1-D FDTD code with PEC boundaries ;7=JU^@D@  
%*********************************************************************** Y>!W&Gtu  
% E#F9<=mA)  
% Program author: Susan C. Hagness tJF~Xv2L!  
% Department of Electrical and Computer Engineering )B5gs%u]  
% University of Wisconsin-Madison 2bG4
,M  
% 1415 Engineering Drive .p*D[o2 9  
% Madison, WI 53706-1691 ``)1`wx$  
% hagness@engr.wisc.edu d`][1rZk  
% %oKc?'L0  
% Copyright 2005 6XCX#4'i%  
% N#!1@!2BN  
% This MATLAB M-file implements a finite-difference time-domain (<~
R[sT|  
% solution of Maxwell's curl equations over a one-dimensional space fy-Z{  
% lattice comprised of uniform grid cells. +6Fdi*:  
% ex!wY  
% To illustrate the algorithm, a 1-GHz sinusoidal wave propagating }eRG$)'  
% in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is _*B~ESC0  
% modeled. The simplified finite difference system for nonpermeable t }C ^E  
% media (discussed in Section 3.6.6 of the text) is implemented. SIVLYi  
% cb&In<q  
% The grid resolution (dx = 1.5 cm) is chosen to provide 20 p:>?  
% samples per wavelength. The Courant factor S=c*dt/dx is set to -oT+;2\2  
% the stability limit (S=1). In 1-D, this is the "magic time step." 8PVs!?Nne  
% The total number of time steps (nmax=240) corresponds to a physical ymVd94L  
% time of 12 ns. `Ta(P30  
% 4O"kOEkKT>  
% The grid is terminated with electric-field components at the far-left 2o}G<7r  
% (i=1) and far-right (i=ie) boundaries. The sinusoidal wave is Z/UVKJm>:  
% launched by an electric-field hard-source condition at i=1 (see 6uE1&-:L  
% Eq. 5.1 in the text). The simplest radiation boundary condition l*MUDT@M8\  
% for plane wave propagation is used to update the electric field ^* v{t?u  
% at i=ie: ]~eWr2uG?  
% q@Yt`$VTN  
% Ez(ie,n+1) = Ez(ie-1,n) i1\ /\^  
% 9nAK6$/  
% To execute this M-file, type "fdtd1D" at the MATLAB prompt. 6i=wAkn_J  
% HD^~4\%  
% This code has been tested in the following Matlab environments: gJ~*rWBK:  
% Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) R6o<p<fTh  
% Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) s31_3?Vdf,  
% Matlab version 7.0.0.19920 R14 (May 6, 2004) :q[n1 O[Ch  
% Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) uB"m!dL  
% Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) mqc Z3lsv  
% dnc!=Z89  
%*********************************************************************** ]=VI"v<X  
{h+E&u[zL  
clear / H/Ne )r  
"/h"Xg>q  
%*********************************************************************** k3h53QTmC  
% Fundamental constants :=7;P)  
%*********************************************************************** q x }fn/:  
/1 %0A  
cc=2.99792458e8; %speed of light in free space TUYl><F5v=  
muz=4.0*pi*1.0e-7; %permeability of free space YH
tI%  
epsz=1.0/(cc*cc*muz); %permittivity of free space o5@P>\u>  
`%I{l  
freq=1.0e+9; %frequency of source excitation vX24W*7  
lambda=cc/freq; %wavelength of source excitation aT1W]i  
omega=2.0*pi*freq; /4BXF4ksi,  
3t6'5{  
%*********************************************************************** `l#$l3v+  
% Grid parameters CP#MNNvgrw  
%*********************************************************************** r>@/XYK&\  
yj9gN}+  
ie=200; %number of Ez samples in grid !+ hgKZ]  
ih=ie-1; %number of Hy samples in grid )L("t  
717m.t,x  
dx=lambda/20.0; %space increment of 1-D lattice T0)y5
  
dt=dx/cc; %time step (S=1.0) [zd-=.:+M[  
omegadt=omega*dt; /s_$CSiB  
qz SI cI  
nmax=round(12.0e-9/dt); %total number of time steps ~i#xjD5  
y+x>{!pw  
x=linspace(dx,ih*dx,ih); +6-!o,(  
CrQ&-!Eh  
%*********************************************************************** 1OeDWEcB  
% Material parameters Z6ex<[`I  
%*********************************************************************** /kVc7LC  
Nn\\}R  
eps=1.0; v
@SrEmg  
sig=5.0e-3; .]l2)O
lLQ  
'J2P3t  
%*********************************************************************** Vs(Zs[  
% Updating coefficients for space region with nonpermeable media -^"?a]B  
%*********************************************************************** -iX!F~qS,  
aJ@qB9(ZBe  
scfact=dt/muz/dx; /Loe y   
NistW+{<  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); N Uml"  
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); 64s;6=  
R7b*(33  
%*********************************************************************** {XW>:EU'N  
% Field arrays DYl{{L8@  
%*********************************************************************** j"=jK^  
\~1+T  
ez(1:ie)=0.0; @C)h;TR  
hy(1:ih)=0.0; 9xp ;$14  
>p:fWQ6  
%*********************************************************************** U2u>A r  
% BEGIN TIME-STEPPING LOOP [=!MS?-G  
%*********************************************************************** e@VRdhb  
^:j:;\;  
for n=1:nmax #=3]bg  
%/C[\wp81  
%*********************************************************************** U<gw<[>f  
% Update electric fields \H12~=p`B  
%*********************************************************************** EZW?(%b>H  
Y>~zt -  
ez(1)=scfact*sin(omegadt*n); h:90K  
@SB+u+mOS  
rbc=ez(ih); >GR
uS\B  
ez(2:ih)=ca*ez(2:ih)+cb*(hy(2:ih)-hy(1:ih-1)); l!'iLq"K(  
ez(ie)=rbc; [!
BH3J!  
Ip-jqN J~  
%*********************************************************************** ,Ou)F;r  
% Update magnetic fields F9hWB17u  
%*********************************************************************** ((q(Q9(F  
go5!zSs  
hy(1:ih)=hy(1:ih)+ez(2:ie)-ez(1:ih); [AwE  
g"f^YEQ_  
%*********************************************************************** __npX_4%S  
% Visualize fields t$|6}BX  
%*********************************************************************** 8~>3&jX  
dj]N59<  
rtime=num2str(round(n*dt/1.0e-9)); @SXgaWr  
gH.^NO5\'  
subplot(2,1,1),plot(x,ez(1:ih)/scfact,'r'),axis([0 3 -1 1]); rP_
)*)  
title(['time = ',rtime,' ns']); ylabel('E_z'); ~`.%n7  
subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]); {`55nwd  
xlabel('x (meters)');ylabel('H_y'); 'M/&bu r  
y9<Fv|Ric  
pause(0.05) C(hg"_W ou  
$ 7!GA9Bn  
%*********************************************************************** N*W.V,6yH  
% END TIME-STEPPING LOOP gQwmYe  
%*********************************************************************** $Yc9><i  
T]`" Xl8  
end u6RHn;b  
x N)Ck76  
我想问的是， @I:&ozy }=  
1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ gV BV@
v!W  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 5Bk  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。 u>Hx#R<*%  
有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把 .. w ?aLWySYT  
mD3#$E!A1  

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

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

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

想看什么样的结果？

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

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

恩，听你的解释，好像我也搞混了。 C5k\RS9  
是看的上面的PEC边界才接着这么理解的。。。 原来下面的编程和PEC没关系？而是吸收边界的处理。 l.gt+
e  
我再回去看看。谢谢了。

我是想看真空中的效果。 IG}`~% Z  
我把把EPS,SIGMA,去掉，并且ca=1,cb=scfact*(dt/epsz),但是出来的图形完全不对。 q>JW$8  
请问问题出在哪呢？ 是作图那边要改么？

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

那个不是相对介质参数的么。 3' Ht
T   
或者改为 eps=1.0; sigma=0;？貌似这样出来的也不对、 aWp9K+4R$/  
还有原来的ca，cb不是也需要修改？ -uho;  
pt"yJtM'P  
那请教怎么改，是真空下的呢？ xB&kxW.;  
我就是想看到，没有衰减的波，遇到边界反射的情况。谢谢了。