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

原程序就是那个一维PEC边界，正弦波在介质中传播的。 ;x:rZV/  
%*********************************************************************** UTf9S>HS  
% 1-D FDTD code with PEC boundaries HcedE3Rg  
%*********************************************************************** H"C[&r  
% 6I!7c^]t  
% Program author: Susan C. Hagness I!>\#K  
% Department of Electrical and Computer Engineering }';D]c  
% University of Wisconsin-Madison W
-
 
% 1415 Engineering Drive HAv{R!*  
% Madison, WI 53706-1691 $2M#qkik-  
% hagness@engr.wisc.edu nK$X[KrV'  
% K-f1{ 0  
% Copyright 2005 FL8g5I  
% om |"S  
% This MATLAB M-file implements a finite-difference time-domain ZR$'u%+g'  
% solution of Maxwell's curl equations over a one-dimensional space Gg~QAsks   
% lattice comprised of uniform grid cells. 5*E#*H  
% rHf&:~   
% To illustrate the algorithm, a 1-GHz sinusoidal wave propagating CBDG./  
% in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is u{FDdR9<  
% modeled. The simplified finite difference system for nonpermeable +<}0|Xl&  
% media (discussed in Section 3.6.6 of the text) is implemented. FJ %  
% p|Q*5
TO  
% The grid resolution (dx = 1.5 cm) is chosen to provide 20 <)\y#N  
% samples per wavelength. The Courant factor S=c*dt/dx is set to KAsS[  
% the stability limit (S=1). In 1-D, this is the "magic time step." {q<03d~9|G  
% The total number of time steps (nmax=240) corresponds to a physical A8q;q2  
% time of 12 ns. <L{(Mj%Z  
% Iw<jT|y)  
% The grid is terminated with electric-field components at the far-left f8SL3+v  
% (i=1) and far-right (i=ie) boundaries. The sinusoidal wave is lip[n;Ir>  
% launched by an electric-field hard-source condition at i=1 (see l3Lyea:  
% Eq. 5.1 in the text). The simplest radiation boundary condition h.!}3\Y  
% for plane wave propagation is used to update the electric field hzI|A~MFB  
% at i=ie: ALEnI@0  
% f>s?4  
% Ez(ie,n+1) = Ez(ie-1,n) S.Z9$k%  
% e#,~,W.H  
% To execute this M-file, type "fdtd1D" at the MATLAB prompt. Fbu5PWhlc  
% PG8^.)]M  
% This code has been tested in the following Matlab environments: =1P6
Vk  
% Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) dB+N\HBY  
% Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) y$3;$ R^  
% Matlab version 7.0.0.19920 R14 (May 6, 2004) Y3h/~bM%
 
% Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) .`7cBsXH  
% Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) K"uNxZ  
% McoK@q;  
%*********************************************************************** {CR5K9  
i 9g>9  
clear R
Jy=pNztm  
>Bs#Xb_B]  
%*********************************************************************** \o\nr!=k  
% Fundamental constants V97,1`  
%***********************************************************************  49d@!  
it>r+%  
cc=2.99792458e8; %speed of light in free space }gkM^*$:%  
muz=4.0*pi*1.0e-7; %permeability of free space \`, [)`  
epsz=1.0/(cc*cc*muz); %permittivity of free space vsL[*OeI  
wBQF~WY  
freq=1.0e+9; %frequency of source excitation &QG6!`fK}3  
lambda=cc/freq; %wavelength of source excitation q %0Cg=  
omega=2.0*pi*freq; {QbvR*gv  
t+?P^Ok  
%*********************************************************************** LTJc,3\,  
% Grid parameters $7QoMV8V  
%*********************************************************************** Q#
(GI2F2#  
RNPbH
.  
ie=200; %number of Ez samples in grid tTN?r8  
ih=ie-1; %number of Hy samples in grid +n
})Y  
+[J/Zw0{  
dx=lambda/20.0; %space increment of 1-D lattice g~BoFc.V2~  
dt=dx/cc; %time step (S=1.0) P/JK$nb  
omegadt=omega*dt; pxSX#S6I  
84i_k  
nmax=round(12.0e-9/dt); %total number of time steps v`V7OD#:j]  
AH4EtZC=W  
x=linspace(dx,ih*dx,ih); JbO ~n )%x  
_SACqamo5s  
%*********************************************************************** |!q$_at  
% Material parameters '!P"xBVAu  
%*********************************************************************** cst}Ibfi  
`VQb-V  
eps=1.0; O)kgBrB  
sig=5.0e-3; f'q 28lVf  
>xA),^ YT  
%*********************************************************************** Z?J:$of*  
% Updating coefficients for space region with nonpermeable media >/<:Q  &  
%*********************************************************************** dr{y0`CCN  
yAL1O94  
scfact=dt/muz/dx; _MWM;f`b  
|wox1Wt|E  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); }X;U|]d  
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); vG^#Sfgtw  
pPVRsXy  
%*********************************************************************** 4yaxl\2  
% Field arrays 5]1leT  
%*********************************************************************** '!Gs>T+  
1!p/6  
ez(1:ie)=0.0; $KjTa#[RX7  
hy(1:ih)=0.0; -ng=l;  
XT,#g-oi  
%*********************************************************************** XYx6V  
% BEGIN TIME-STEPPING LOOP 4$,,Ppn  
%*********************************************************************** ha;l(U>  
AGYm';z3  
for n=1:nmax Ufo>|A6;$  
BpO9As 1um  
%*********************************************************************** *m+5Pr`7  
% Update electric fields K?#]("De6  
%*********************************************************************** XE}H3/2  
Ip}Vb6}  
ez(1)=scfact*sin(omegadt*n); 4z:#I;  
rZ_>`}O2  
rbc=ez(ih); gQ~5M'#  
ez(2:ih)=ca*ez(2:ih)+cb*(hy(2:ih)-hy(1:ih-1)); IfDx@?OB  
ez(ie)=rbc; +VEU:1Gt  
">dq0gD  
%*********************************************************************** RV-hIdAU  
% Update magnetic fields OlxX.wP  
%*********************************************************************** -Uo?WXP]B'  
*jzLFuWIG  
hy(1:ih)=hy(1:ih)+ez(2:ie)-ez(1:ih); ovf/;Q/}  
LF*Q!  
%*********************************************************************** pz_e=xr  
% Visualize fields ^_p%Yv  
%*********************************************************************** A1cb"N^  
ly4Qg\l  
rtime=num2str(round(n*dt/1.0e-9)); Jf:,y~mV  
-%IcYz
yA  
subplot(2,1,1),plot(x,ez(1:ih)/scfact,'r'),axis([0 3 -1 1]); yy2Ie  
title(['time = ',rtime,' ns']); ylabel('E_z'); FM^9
}*  
subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]); 
Gie@JX  
xlabel('x (meters)');ylabel('H_y'); Y9r3XhVI  
y[0`hSQ)~  
pause(0.05) lm'Zy"~::  
=l
r)gj  
%*********************************************************************** IGj`_a  
% END TIME-STEPPING LOOP ~@I@}n  
%*********************************************************************** S+x_c4 T  
sCH)gr@gJ^  
end 4ax|Vb)D  
Uhh l3%p  
我想问的是， H5wb_yBQ+  
1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ `?s.\Dh  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 7CvD'QW /  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。 ~el-*=<m  
有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把 .. =N.!k Vkl  
v:ER4  

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

1、PEC就是理想导体边界，是程序里是不是rbc=ez(ih); ez(ie)=rbc;这2句。代替了“i=ie处， Ez(ie,n+1) = Ez(ie-1,n)”这个？ &7][@v  
2、理想导体边界的话，那不是应该是全反射？但是，看到的波形只是衰减的传播。。并没有看到反射啊。 [$:@X V(  
这两句是吸收电磁波的，所以加了以后应该没有反射波了。 dpn3 (  
tR O IBq|  
3、我把rbc=ez(ih); ez(ie)=rbc，这2句边界条件去了，看到的反而有反射。有点迷糊。在我的理解，应该是有这个条件的话，应该看到的是有反射，把那个去了，才看到的是没有反射，继续传播啊？ csvOg[  
对的，去掉以后，边界就相当于电场为0，等同于金属边界了。

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

想看什么样的结果？

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

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

恩，听你的解释，好像我也搞混了。 KlSY^(kHR  
是看的上面的PEC边界才接着这么理解的。。。 原来下面的编程和PEC没关系？而是吸收边界的处理。 'DB({s  
我再回去看看。谢谢了。

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

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

那个不是相对介质参数的么。 n[-!Jp[  
或者改为 eps=1.0; sigma=0;？貌似这样出来的也不对、 d8I:F9  
还有原来的ca，cb不是也需要修改？ I(S6DkU  
}e@-
[RJ!  
那请教怎么改，是真空下的呢？ |`/uS;O  
我就是想看到，没有衰减的波，遇到边界反射的情况。谢谢了。