[求助]请高手进来指点一下，小人用matlab编的FDTD的PML吸收边界总是出问题

那位能人能能帮看下以下程序将不胜感激，只要帮我看看左上角和左边界的PML吸收条件就可以了，谢谢各位大侠了。 +*OAClt+]  
function [Ez,haoshi,Hx,Hy]=PML(xunhuan)   %TE model Le83[E*i  
t0=clock; 0 Rb3|te  
d=9; %C%3c4+Oh  
N=121+2*d; D)MFii1J~  
Ez=zeros(N,N); i59}6u_f  
Hx=zeros(N,N); A":=-$)  
Hy=zeros(N,N); Q``1^E'  
Ezt=zeros(N,N); )]n>.ZmLCB  
Hxt=zeros(N,N); N_G&nw  
Hyt=zeros(N,N); G!%m~+",  
"=2\kZ  
Ezxt=zeros(N,N); j6s j2D  
Ezyt=zeros(N,N); LYAG
pcG  
Ezx=zeros(N,N); x,: k/]  
Ezy=zeros(N,N); 2fdN@iruB  
DXJw)%G w  
ox=zeros(N,N); Dx\~#$S!=  
oy=zeros(N,N); Zzlt^#KLx  
omx=zeros(N,N); 5H'Iul<Os  
omy=zeros(N,N); ZQ"dAR/y  
@l3&vt2=J  
*vQ 6LF;y  
9z?c0W5x  
Qfp4}a=  
%--------------------------------常数设定----------------------------------- FM%WMyb[  
lamda=7e-6;                                                                %波长(m) 7!qeIz  
time=xunhuan;                                                                 %时间数 a|TUH+|  
dx=lamda/14;                                                                %横轴间距 =nHkFi@D=t  
dy=lamda/14;                                                                %纵轴间距 E2l"e?AN~  
c=3e+8;                                                                    %光速(m) #
@nPB.  
w=2*pi*c/lamda;                                                            %角频率             eP (*.  
T=2*pi/w;                                                                  %周期 }S"gZ6  
dt=T/40;                                                                   %时间间距 L@RnLaoQ  
e0=8.85e-12; aGWO3Nk  
u0=4e-7*pi; 6l,6k~Z9  
m=4; -qpvV
LR,  
JQLQS  
omax=-(m+1)*log(1e-2)/2/sqrt(u0/e0)/d; 46M=R-7=  
ommax=u0*omax/e0; e O~p"d-|  
/1TK+E$  
j7d^ga-`  
ZKW1HL ]m  
3 85qQppz  
;\"5)S  
I=10:130; 
_I)TO_L;  
J=10:130; 3N]ushMO  
8-gl$h  
I1=2:9; 5xT, O  
I2=131:138; &pY$\  
i=5!taxu}E  
J1=2:9; "[/W+&z[~  
J2=131:138; ,or;8aYc#  
T6SYXQd>.  
@Y,t]  
for p=1:9 ^?|4<Rm  
    ox(:,10-p)=omax*(p/d)^m; =f@
71D1  
    omx(:,10-p)=ommax*(p/d)^m; f2tCB1[D+  
    ox(:,130+p)=omax*(p/d)^m; J_a2DM6d  
    omx(:,130+p)=ommax*(p/d)^m; zM8 jjB  
     LQqba4$  
     |5(CzXR]  
    oy(10-p,:)=omax*(p/d)^m; nVlZ_72d  
    omy(10-p,:)=ommax*(p/d)^m; };rEN`L  
    oy(130+p,:)=omax*(p/d)^m; l_2Xao$  
    omy(130+p,:)=ommax*(p/d)^m;     %-YWn`yEm  
end wBlE!Pm  
K}q5,P(  
K1=dt/dy/e0; "z6p=B"?3  
K2=dt/dy/u0; C
m5L99Y  
{%6 '|<`[  
Vh01y f  
+dCR$<e9
r  
Z$i?p;HnW  
@@a#DjE%/  
[6\O <-?  
for n=1:time -^np"Jk  
%左上角 I1, J1 Na]ITCVR  
% % ;"f9"  
Hx(I1,J1)=exp(-omy(I1,J1).*dt/u0).*Hxt(I1,J1)-(1-exp(-omy(I1,J1).*dt/u0))./(dy.*omy(I1,J1)).*(Ezt(I1,J1+1)-Ezt(I1,J1)); d~.hp  
Hy(I1,J1)=exp(-omx(I1,J1).*dt/u0).*Hyt(I1,J1)+(1-exp(-omy(I1,J1).*dt/u0))./(dx.*omx(I1,J1)).*(Ezt(I1+1,J1)-Ezt(I1,J1)); pVl7]_=m  
Ezx(I1,J1)=exp(-ox(I1,J1).*dt/e0).*Ezxt(I1,J1)+(1-exp(-ox(I1,J1).*dt/e0))./(dx.*ox(I1,J1)).*(Hy(I1,J1)-Hy(I1-1,J1)); Gqq<-drR  
Ezy(I1,J1)=exp(-oy(I1,J1).*dt/e0).*Ezyt(I1,J1)-(1-exp(-oy(I1,J1).*dt/e0))./(dy.*oy(I1,J1)).*(Hx(I1,J1)-Hx(I1,J1-1)); T&1-gswr:  
% io"NqR#"v  
Ezyt(I1,J1)=Ezy(I1,J1); Ezxt(I1,J1)=Ezx(I1,J1); Ez(I1,J1)=Ezx(I1,J1)+Ezy(I1,J1); Hxt(I1,J1)=Hx(I1,J1); Hyt(I1,J1)=Hy(I1,J1); a^G>|+8  
Ezt(I1,J1)=Ez(I1,J1); DZ`,QWuA  
']Czn._  
% dM,{:eID  
% 上边界 I1, J +U'n|>t9  
c.Izm+9k  
Hx(I1,J)=exp(-omy(I1,J).*dt/u0).*Hxt(I1,J)-(1-exp(-omy(I1,J).*dt/u0))./(dy.*omy(I1,J)).*(Ezt(I1,J+1)-Ezt(I1,J)); A?-t`J  
Hy(I1,J)=Hyt(I1,J)+K2.*(Ezt(I1+1,J)-Ezt(I1,J)); I+Y Z+  
Ezx(I1,J)=Ezxt(I1,J)+K1.*(Hy(I1,J)-Hy(I1-1,J)); }ub>4N[  
Ezy(I1,J)=exp(-oy(I1,J).*dt/e0).*Ezyt(I1,J)-(1-exp(-oy(I1,J).*dt/e0))./(dy.*oy(I1,J)).*(Hx(I1,J)-Hx(I1,J-1)); oGXcu?ft  
s!W{ru  
C(sz/x?11  
1C=42ZZ&2  
Ezyt(I1,J)=Ezy(I1,J); Ezxt(I1,J)=Ezx(I1,J); Ez(I1,J)=Ezx(I1,J)+Ezy(I1,J); Hxt(I1,J)=Hx(I1,J); Hyt(I1,J)=Hy(I1,J); z$Z%us>io  
Ezt(I1,J)=Ez(I1,J); JFu.o8[Q  
R {-M%n4w  
% %右上角 I1, J2 +pUYFDwFx  
% dr8Q>(ZY  
% * W"Pv,:  
% Hx(I1,J2)=exp(-omy(I1,J2).*dt/u0).*Hxt(I1,J2)-(1-exp(-omy(I1,J2).*dt/u0))./(dy.*omy(I1,J2)).*(Ezt(I1,J2+1)-Ezt(I1,J2)); ! qtj1.w  
% Hy(I1,J2)=exp(-omx(I1,J2).*dt/u0).*Hyt(I1,J2)+(1-exp(-omy(I1,J2).*dt/u0))./(dx.*omx(I1,J2)).*(Ezt(I1+1,J2)-Ezt(I1,J2)); mE+=H]`.p  
% Ezx(I1,J2)=exp(ox(I1,J2).*dt/e0).*Ezxt(I1,J2)+(1-exp(ox(I1,J2).*dt/e0))./(dx.*ox(I1,J2)).*(Hy(I1,J2)-Hy(I1-1,J2)); Yjy%
MR  
% Ezy(I1,J2)=exp(-oy(I1,J2).*dt/e0).*Ezyt(I1,J2)-(1-exp(-oy(I1,J2).*dt/e0))./(dy.*oy(I1,J2)).*(Hx(I1,J2)-Hx(I1,J2-1)); W`[7|8(6!  
% % xt@v"P2Ok  
% Ezyt(I1,J2)=Ezy(I1,J2); Ezxt(I1,J2)=Ezx(I1,J2); Ez(I1,J2)=Ezx(I1,J2)+Ezy(I1,J2); Hxt(I1,J2)=Hx(I1,J2); Hyt(I1,J2)=Hy(I1,J2); E-X02A  
% Ezt(I1,J2)=Ez(I1,J2); 7lOAu]Zx  
% % %左边界 I,J1 _x-2tnIxXv  
Hx(I,J1)=Hxt(I,J1)-K2.*(Ezt(I,J1+1)-Ezt(I,J1)); ,LOx!  
Hy(I,J1)=exp(-omx(I,J1).*dt/u0).*Hyt(I,J1)+(1-exp(-omy(I,J1).*dt/u0))./(dx.*omx(I,J1)).*(Ezt(I+1,J1)-Ezt(I,J1)); \HMuVg'Q  
Ezx(I,J1)=exp(-ox(I,J1).*dt/e0).*Ezxt(I,J1)+(1-exp(-ox(I,J1).*dt/e0))./(dx.*ox(I,J1)).*(Hy(I,J1)-Hy(I-1,J1)); pcd?6jh8  
Ezy(I,J1)=Ezyt(I,J1)-K1.*(Hx(I,J1)-Hx(I,J1-1)); yqJ>Z%)hf  
_4{3^QZq5  
Ezyt(I,J1)=Ezy(I,J1); Ezxt(I,J1)=Ezx(I,J1); Ez(I,J1)=Ezx(I,J1)+Ezy(I,J1); Hxt(I,J1)=Hx(I,J1); Hyt(I,J1)=Hy(I,J1); i*xVD`x~  
Ezt(I,J1)=Ez(I,J1); C
9Cl$yZ  
%中心循环区域 x wfdJ(&  
Ez(I,J)=Ezt(I,J)+dt/e0*(-((Hxt(I,J)-Hxt(I,J-1))/dy-(Hyt(I,J)-Hyt(I-1,J))/dx)); K-(C5 "j_  
Ez(70,70)=sin(w*(n+1/2)*dt)*kaiguan((n+1/2)*dt,T);                            %光源设定 5a5JOl$8  
Hy(I,J)=Hyt(I,J)+dt/u0*((Ez(I+1,J)-Ez(I,J))/dx); lBG5~<NT  
Hx(I,J)=Hxt(I,J)+dt/u0*(-(Ez(I,J+1)-Ez(I,J))/dy); @Us#c 7/  
Ezt(I,J)=Ez(I,J); Hxt(I,J)=Hx(I,J); Hyt(I,J)=Hy(I,J);                                                 %上次时间场的分布 # b3 14  
% I_c?Ky8J_|  
% %右边界 I, J2
{6}$XLV3l  
Hx(I,J2)=Hxt(I,J2)-K2.*(Ezt(I,J2+1)-Ezt(I,J2)); RR`\q>|  
Hy(I,J2)=exp(-omx(I,J2).*dt/u0).*Hyt(I,J2)+(1-exp(-omy(I,J2).*dt/u0))./(dx.*omx(I,J2)).*(Ezt(I+1,J2)-Ezt(I,J2)); %Zeb#//Jz  
Ezx(I,J2)=exp(-ox(I,J2).*dt/e0).*Ezxt(I,J2)+(1-exp(-ox(I,J2).*dt/e0))./(dx.*ox(I,J2)).*(Hy(I,J2)-Hy(I-1,J2)); M6[O>z  
Ezy(I,J2)=Ezyt(I,J2)-K1.*(Hx(I,J2)-Hx(I,J2-1)); x {Rj2~KC  
qP/McH?  
Ezyt(I,J2)=Ezy(I,J2); Ezxt(I,J2)=Ezx(I,J2); Ez(I,J2)=Ezx(I,J2)+Ezy(I,J2); Hxt(I,J2)=Hx(I,J2); Hyt(I,J2)=Hy(I,J2); 3;//o<  
Ezt(I,J2)=Ez(I,J2); AAi4} 8+\  
% %左下角 I2, J1 %@I= $8j  
% Hx(I2,J1)=exp(-omy(I2,J1).*dt/u0).*Hxt(I2,J1)-(1-exp(-omy(I2,J1).*dt/u0))./(dy.*omy(I2,J1)).*(Ezt(I2,J1+1)-Ezt(I2,J1)); )Zvn{  
% Hy(I2,J1)=exp(-omx(I2,J1).*dt/u0).*Hyt(I2+1,J1)-(1-exp(-omy(I2,J1).*dt/u0))./(dx.*omx(I2,J1)).*(Ezt(I2,J1)-Ezt(I2+1,J1));                          %%%% *P12d  
% Ezx(I2,J1)=exp(-ox(I2,J1).*dt/e0).*Ezxt(I2,J1)-(1-exp(-ox(I2,J1).*dt/e0))./(dx.*ox(I2,J1)).*(Hy(I2-1,J1)-Hy(I2,J1)); 1<]?@[l<  
% Ezy(I2,J1)=exp(-oy(I2,J1).*dt/e0).*Ezyt(I2,J1)-(1-exp(-oy(I2,J1).*dt/e0))./(dy.*oy(I2,J1)).*(Hx(I2,J1)-Hx(I2,J1-1)); I'J-)D`  
% +nYF9z2  
% Ezyt(I2,J1)=Ezy(I2,J1); Ezxt(I2,J1)=Ezx(I2,J1); Ez(I2,J1)=Ezx(I2,J1)+Ezy(I2,J1); Hxt(I2,J1)=Hx(I2,J1); Hyt(I2,J1)=Hy(I2,J1); eH;{Ln  
% Ezt(I2,J1)=Ez(I2,J1); vA?3kfL|#  
% %下边界 I2, J 6m9\0)R  
Hx(I2,J)=exp(-omy(I2,J).*dt/u0).*Hxt(I2,J)-(1-exp(-omy(I2,J).*dt/u0))./(dy.*omy(I2,J)).*(Ezt(I2,J+1)-Ezt(I2,J)); .%\R L/  
Hy(I2,J)=Hyt(I2,J)+K2.*(Ezt(I2+1,J)-Ezt(I2,J)); ![[:Z  
Ezx(I2,J)=Ezxt(I2,J)+K1.*(Hy(I2,J)-Hy(I2-1,J)); Z'wGZ(  
Ezy(I2,J)=exp(-oy(I2,J).*dt/e0).*Ezyt(I2,J)-(1-exp(-oy(I2,J).*dt/e0))./(dy.*oy(I2,J)).*(Hx(I2,J)-Hx(I2,J-1)); #E/|WT  
lo7>
$`Q  
Ezyt(I2,J)=Ezy(I2,J); Ezxt(I2,J)=Ezx(I2,J); Ez(I2,J)=Ezx(I2,J)+Ezy(I2,J); Hxt(I2,J)=Hx(I2,J); Hyt(I2,J)=Hy(I2,J); !4"$O@U4  
Ezt(I2,J)=Ez(I2,J); BgsU:eKe  
% %右下角  I2，J2 f1\mE~#}  
% Hx(I2,J2)=exp(-omy(I2,J2).*dt/u0).*Hxt(I2,J2+1)-(1-exp(-omy(I2,J2).*dt/u0))./(dy.*omy(I2,J2)).*(Ezt(I2,J2+1)-Ezt(I2,J2)); 9:!V":8q  
% Hy(I2,J2)=exp(-omx(I2,J2).*dt/u0).*Hyt(I2+1,J2)-(1-exp(-omy(I2,J2).*dt/u0))./(dx.*omx(I2,J2)).*(Ezt(I2,J2)-Ezt(I2+1,J2));                          %%%% JWQd6JQ_~V  
% Ezx(I2,J2)=exp(-ox(I2,J2).*dt/e0).*Ezxt(I2,J2)-(1-exp(-ox(I2,J2).*dt/e0))./(dx.*ox(I2,J2)).*(Hy(I2-1,J2)-Hy(I2,J2)); w\JTMS$  
% Ezy(I2,J2)=exp(-oy(I2,J2).*dt/e0).*Ezyt(I2,J2)-(1-exp(-oy(I2,J2).*dt/e0))./(dy.*oy(I2,J2)).*(Hx(I .. e9F+R@8  
|UQGZ  

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

没人理呀…………自己来顶

二维的怎么不用mur吸收边界呢，效果也挺好的！

mur用过了，想试下PML

http://bbs.rfeda.cn/read-htm-tid-16687-fpage-0-toread--page-2.html sm/aL^4  
这里面有个源码PML+2D FDTD。参考一下吧。