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

那位能人能能帮看下以下程序将不胜感激，只要帮我看看左上角和左边界的PML吸收条件就可以了，谢谢各位大侠了。 1p/_U?H:|  
function [Ez,haoshi,Hx,Hy]=PML(xunhuan)   %TE model $*XTX?,'  
t0=clock;
^^uY)AL  
d=9; [<^'}-SJ  
N=121+2*d; 1z,P"?Q  
Ez=zeros(N,N); 3-;<G  
Hx=zeros(N,N); -c0*  
Hy=zeros(N,N); 3S>rc0]6  
Ezt=zeros(N,N); QDK }e:4q  
Hxt=zeros(N,N); 9}K K]m6u}  
Hyt=zeros(N,N); GX.a!XQ@!  
.hf%L1N%F  
Ezxt=zeros(N,N); n sN n>{  
Ezyt=zeros(N,N); y^X]q[-?  
Ezx=zeros(N,N); f@Ve,i  
Ezy=zeros(N,N); ~S :8M<aB  
K-.%1d@$y  
ox=zeros(N,N); #QOb[9(Tu(  
oy=zeros(N,N); 8 f~M
6  
omx=zeros(N,N); $_a/!)bP  
omy=zeros(N,N); ]w-W  
??hKsjNAm0  
$S?xB$  
mH<|.7~0  
IK4(r /  
%--------------------------------常数设定----------------------------------- hf)RPG&  
lamda=7e-6;                                                                %波长(m) ]E.FBGT  
time=xunhuan;                                                                 %时间数 V^;lg[:  
dx=lamda/14;                                                                %横轴间距 NDe FY  
dy=lamda/14;                                                                %纵轴间距 *r~6R
  
c=3e+8;                                                                    %光速(m) 7oL:C  
w=2*pi*c/lamda;                                                            %角频率             H-m).^  
T=2*pi/w;                                                                  %周期 j\BtaC  
dt=T/40;                                                                   %时间间距 ,&O&h2=  
e0=8.85e-12; i3(5 '  
u0=4e-7*pi; HyQ(9cn|  
m=4;
$b_~  
6|6O| <o  
omax=-(m+1)*log(1e-2)/2/sqrt(u0/e0)/d; G%jV}
7h  
ommax=u0*omax/e0; eXLdb-  
]P^3uXi  
}LWrtmc  
ja{x}n*5  
=zp{ ^mC  
cqb6]  
I=10:130; m+pK,D~{"  
J=10:130; @@!]Raj=  
.o<9[d"  
I1=2:9; o/RGzPR  
I2=131:138; 2of+KI:
 
ay{]Vqi9  
J1=2:9; _39VL  
J2=131:138; Q"LlBp>t|#  
s9u7zqCF  
sG|,#XQ  
for p=1:9 Z#;\Rb.x7  
    ox(:,10-p)=omax*(p/d)^m; Ym-mfWo^#  
    omx(:,10-p)=ommax*(p/d)^m; !.q#X^@>L  
    ox(:,130+p)=omax*(p/d)^m; S =sL:FC  
    omx(:,130+p)=ommax*(p/d)^m; p .~5k  
     i MS4<`  
     0J5$ Yw1'F  
    oy(10-p,:)=omax*(p/d)^m; S~g"  
    omy(10-p,:)=ommax*(p/d)^m; "zIQ(|TL?d  
    oy(130+p,:)=omax*(p/d)^m; bg|=)sw4  
    omy(130+p,:)=ommax*(p/d)^m;     !0X"^VB  
end WUx2CK2N  
Kt"4<'  
K1=dt/dy/e0; UG
]5Dxk  
K2=dt/dy/u0; _Mh..#)`[  
^(qR({cX  
r6:nYyF$)v  
5RSP.Vyx{  
bGj<Dojl  
pqbKPpG  
jlD3SF~2  
for n=1:time <V8=*n"mR  
%左上角 I1, J1 7Z81+I|&8  
% % MLDAr dvK  
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)); 7B)@ aUj$  
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)); ;O .;i,#Z  
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)); THwq~c'  
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)); M?ElD1#Z  
% fL&e^Q  
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); Z= pvoTY  
Ezt(I1,J1)=Ez(I1,J1); bh5C  
BJZGQrsz  
% 2m&?t_W  
% 上边界 I1, J :|ytw=3>  
K}LF ${bS  
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)); iA"H*0  
Hy(I1,J)=Hyt(I1,J)+K2.*(Ezt(I1+1,J)-Ezt(I1,J)); 0}Qd  
Ezx(I1,J)=Ezxt(I1,J)+K1.*(Hy(I1,J)-Hy(I1-1,J));  t|:XSJ9  
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)); 'GZ,  
E3_ 5~>  
$A:?o?"7}  
bde6 ;=oM  
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); 5XNFu C9E  
Ezt(I1,J)=Ez(I1,J); ugW.nf*O  
QP6a,^];  
% %右上角 I1, J2 c 8|&Q  
% aQ1n1OBr  
% {\k:?w4  
% 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)); O;#0Yg  
% 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)); R"71)ob4  
% 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)); <$nMqUu0  
% 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)); OZl0I#@A  
% % lYrW"(2  
% 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); '&y+,2?;Y[  
% Ezt(I1,J2)=Ez(I1,J2); M;0\fUh;  
% % %左边界 I,J1 lIatM@gU  
Hx(I,J1)=Hxt(I,J1)-K2.*(Ezt(I,J1+1)-Ezt(I,J1)); l_&T)Ei  
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)); 7<F{a"5P  
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));  !h*F58  
Ezy(I,J1)=Ezyt(I,J1)-K1.*(Hx(I,J1)-Hx(I,J1-1)); E{B40E~4  
c/G^}d%  
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); 4e|(= W`  
Ezt(I,J1)=Ez(I,J1); QnH~'
 k  
%中心循环区域 z-kB!~r  
Ez(I,J)=Ezt(I,J)+dt/e0*(-((Hxt(I,J)-Hxt(I,J-1))/dy-(Hyt(I,J)-Hyt(I-1,J))/dx)); SYv5{bff =  
Ez(70,70)=sin(w*(n+1/2)*dt)*kaiguan((n+1/2)*dt,T);                            %光源设定 X5P1wxk'  
Hy(I,J)=Hyt(I,J)+dt/u0*((Ez(I+1,J)-Ez(I,J))/dx); m8v=pab e  
Hx(I,J)=Hxt(I,J)+dt/u0*(-(Ez(I,J+1)-Ez(I,J))/dy); 3.04Toq!  
Ezt(I,J)=Ez(I,J); Hxt(I,J)=Hx(I,J); Hyt(I,J)=Hy(I,J);                                                 %上次时间场的分布 O~F8lQ  
% =I)Ex)  
% %右边界 I, J2 VZU@G)rd  
Hx(I,J2)=Hxt(I,J2)-K2.*(Ezt(I,J2+1)-Ezt(I,J2)); 9q<?
xO  
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)); *3y:Wv T>  
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)); ur/:aI  
Ezy(I,J2)=Ezyt(I,J2)-K1.*(Hx(I,J2)-Hx(I,J2-1)); h{VGhkU9f  
$
K~ t'wr  
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); $Mqw)X&q  
Ezt(I,J2)=Ez(I,J2); @0%^\Qf2  
% %左下角 I2, J1 S#Pni}JD  
% 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)); ]~m2#g%  
% 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));                          %%%% _ 3jY,*  
% 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)); F[oTc^dr  
% 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)); )G$0:-J-  
% tp+H]H3  
% 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); ,09d"7`X  
% Ezt(I2,J1)=Ez(I2,J1); a#P
{[  
% %下边界 I2, J sHMZ'9b  
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<8!lQ4s  
Hy(I2,J)=Hyt(I2,J)+K2.*(Ezt(I2+1,J)-Ezt(I2,J)); +q~dS.  
Ezx(I2,J)=Ezxt(I2,J)+K1.*(Hy(I2,J)-Hy(I2-1,J)); {(`xA,El  
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)); AkV8}>G?#A  
;e^`r;]  
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); v6KF0mqA&  
Ezt(I2,J)=Ez(I2,J); (wEaw|Zx  
% %右下角  I2，J2 ljO t~@Ea  
% 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)); 9iOTT%pq  
% 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));                          %%%% 'uF"O"*  
% 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)); 'Y-Y By :  
% 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 .. :]IYw!_-p  
P:HmT   

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

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

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

mur用过了，想试下PML

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