[已解决]Mur二阶边界条件求助

前面编了一阶的,能很好的运行,现在想编一个二阶的,于是换了个公式,但是当波进入边界后,马上不正常了(出现超过极限的反射),下面就是使用的二阶Mur的边界条件(按葛书编的),望大家指教. Xj+_
"0 #  
% 3.二阶Mur边界条件 K
RlJKd{  
k0=8.854e-12;%真空介电常数 .S|T{DMQ[  
h1=(c0*dt-ds)/(c0*dt+ds); r=3`Eb"t  
h2=c0^2*k0*dt/2/(c0*dt+ds); @mZK[*Ak<*  
7#+Ih-&EQ  
for i=1 Y>aVnixx<  
    for j=2:JE-1 r4[=pfe25  
        Hz(i,j)=Hz(i+1,j)+h1*(Hz(i+1,j)-Hz(i,j))+h2*(Ex(i,j)-Ex(i,j-1)+Ex(i+1,j)-Ex(i+1,j-1));%左侧边界 aNKw.S>  
    end sx azl]
 
end jt}oq%Bf  
for i=IE ]\K?%z
  
    for j=2:JE-1 !I1p`_(_7  
        Hz(i,j)=Hz(i-1,j)+h1*(Hz(i-1,j)-Hz(i,j))+h2*(Ex(i,j)-Ex(i,j-1)+Ex(i-1,j)-Ex(i-1,j-1));%右侧边界 aWimg6q  
    end X\!q8KEpR&  
end Pq<43:*?  
for j=1 pSC{0Y$g  
    for i=2:IE-1 ^PFiO 12  
        Hz(i,j)=Hz(i,j+1)+h1*(Hz(i,j+1)-Hz(i,j))+h2*(Ey(i,j)-Ey(i-1,j)+Ey(i,j+1)-Ey(i-1,j+1));%下侧边界 M:OZW
YQ  
    end @F(e
r  
end gQik>gFr  
for j=JE t)8crX}P  
    for i=2:IE-1 vcy1itY  
        Hz(i,j)=Hz(i,j-1)+h1*(Hz(i,j-1)-Hz(i,j))+h2*(Ey(i,j)-Ey(i-1,j)+Ey(i,j-1)-Ey(i-1,j-1));%上侧边界 KL?<lp"  
 & .. i#YDdz  
"k+ :!D  

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

好像公式不是这样的吧。

找本葛德彪老师那本书后面有个2维2阶Mur的例子，fortran编的，和matlab语法差不多

能不能给我一下matlab的二阶Mur边界条件,我手头上现在没书.

Mur2阶的公式有4个系数，公式有问题

% 本程序实现2维TM波FDTD仿真 N06O.bji  
% 此程序用PML设置吸收边界条件 PM!t"[@&  
% FDTD_2D_PML MWh+h7k'  
% 仅含有Ez,Hx,Hy分量 u#k,G`  
=[os<+  
clear; arVf"3a  
clc; ~j#6 goKn  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M|e n>P  
%% 1.初始化 w]o5L  
T=200; % 迭代次数 y7EX&  
IE=100; % yc=#Jn?S  
JE=100; !&p:=}s  
npml=8; % PML的网格数量 e7qMt[.  
{0WIDD  
c0=3*10^8; % 波速 <JM%Kn )  
f=1.5*10^(9); % 频率 vqz#V=J{  
lambda=c0/f; % 波长 ~{Rt4o _W  
L$c%u
 
wl=10; %数值波长 -G#@BtB2+  
dx=lambda/wl; iiB)/~!O  
dy=lambda/wl; B\>}X_\4  
pi=3.14159; ]G~N+\8]U  
*Ne2l`!1m  
dt=dx/(2*c0); % 时间间隔 ikG9l&n  
JD`;,Md  
epsz=1/(4*pi*9*10^9); % 真空介电常数 }.V0SM6  
epsilon=1; % 相对介电常数 _XNR um4  
sigma=0; % 电导率 `+Z#*lj|@  
qS}RFM5|  
spread=6; % 脉冲宽度 5>[s
Cl-  
t0=20; % 脉冲高度  / !  
RdvTtXg  
ic=IE/2; % 源的X位置 (`.qG &6p  
jc=JE/2; % 源的Y位置 gdFoTcHgO|  
%初始化 bTy)0ta>AF  
for i=1E+1; He-Ja  
for j=1:JE+1; r9a!,^}F  
dz(i,j)=0; % z方向电荷密度 s1 ^mk]  
ez(i,j)=0; % z方向电场 Hs-.83V  
hx(i,j)=0; % x方向磁场 JX$NEq(  
hy(i,j)=0; % y方向磁场 (i0"hi  
ihx(i,j)=0;% as|c`4r\O  
ihy(i,j)=0; u+'@>%7  
iz(i,j)=0; % z方向求和参量，频域卷积转化为时域求和 B "F`OS[  
end; )/T$H|  
end; Yn]yd1  
^npS==Y]!.  
for i=2E; % ^'53]b:  
for j=2:JE; $0S#d@v}  
ga(i,j)=1; ) a\DS yr  
end; `e'o~oSu  
end; 2n;;Tso"  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  [Ro0eH  
%PML参数的设置 '7LJuMp$#  
for i=1E; q '{<c3&  
gi2(i)=1; Fif^V  
gi3(i)=1; <bXWkj  
fi1(i)=0; 2G}7R5``9  
fi2(i)=1.0; A.C278^O8  
fi3(i)=1.0; ;E? hz  
end = *;Xc-_  
v\9,j  
for j=1:JE; Fp6[W5>(-  
gj2(j)=1; 4OZ5hH h  
gj3(j)=1; ,#jhKnk2e  
fj1(j)=0; HZ<f(  
fj2(j)=1; [(hvK{)  
fj3(j)=1; XvkI+c  
end (6#yw`\  
Ed0>
R<jR9  
for i=1:npml+1; %设置PML层中的参数 R ms
01m>Y  
xnum=npml+1-i;  1C,C)  
xn=0.33*(xnum/npml)^3; !0csNg!  
gi2(i)=1.0/(1+xn); LVFsd6:h  
gi2(IE-1-i)=1/(1+xn); a.&#dxgW[  
gi3(i)=(1-xn)/(1+xn); qwNKRqT  
gi3(IE-1-i)=(1-xn)/(1+xn); E<#4G9O<  
xn=0.25*((xnum-0.5)/npml)^3; Xt,,AGm}  
fi1(i)=xn; pg}+lYGP  
fi1(IE-2-i)=xn; +V+*7s%fL  
fi2(i)=1.0/(1+xn); (06Vcqg  
fi2(IE-2-i)=1/(1+xn); ~eV!!38 J  
fi3(i)=(1-xn)/(1+xn); 1QF*e'  
fi3(IE-2-i)=(1-xn)/(1+xn); 6n6VEwYj  
end &;h~JS=  
~UJu @M  
for i=1:npml+1; b~Pxgfu"  
xnum=npml+1-i; 90h1e7Zc
C  
xn=0.33*(xnum/npml)^3; &Wz`>qYL*  
gj2(i)=1.0/(1+xn); &kjwIg{  
gj2(JE-1-i)=1/(1+xn); +x9"#0|k;  
gj3(i)=(1-xn)/(1+xn); Sd<@X@iU8D  
gj3(JE-1-i)=(1-xn)/(1+xn); 9<(K6Q  
xn=0.25*((xnum-0.5)/npml)^3; HMS9y%zl/  
fj1(i)=xn; @+\S!o3m  
fj1(JE-2-i)=xn; H,XLb.  
fj2(i)=1.0/(1+xn); 2hwXWTSu  
fj2(JE-2-i)=1/(1+xn); ci_v7Jnwo  
fj3(i)=(1-xn)/(1+xn); 3+
8"  
fj3(JE-2-i)=(1-xn)/(1+xn);  kulQR>u  
end 'f?&EsIV?  
HLoQ}oK|K  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FH`'1iVH  
%% 2.迭代求解电场和磁场 Y9@dZw%2  
for t=1:T; R|)2Dg  
for i=2E; % 为了使每个电场周围都有磁场进行数组下标处理 B`?N0t%X  
for j=2:JE; 78a
-3){  
dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1)); A?;8%00  
end; /@:up+$  
end; % 电场循环结束 Msa6yD#  
C/CfjRzd  
pulse=sin(2*pi*f*t*dt); % 正弦波源 gR-Qj  
dz(ic,jc)=dz(ic,jc)+pulse; % 软源 X#kjt)W  
~ S?-{
X+  
for i=1E; % 为了使每个电场周围都有磁场进行数组下标处理 h\u0{!@}  
for j=1:JE;
qzHqj;  
ez(i,j)=ga(i,j)* dz(i,j); %反映煤质的情况都是放到这里的 v<7Gln  
% iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ; +{C9uY)$vf  
end; ;&W;  
end; % 电荷密度循环结束 T$8@2[  
^0s\/qyqm  
for j=1:JE; QLB1:O>  
ez(1,j)=0; g<rKV+$6  
ez(IE,j)=0; :B*vkwT  
end ^QXw[th!d  
zOiY0`=  
for i=1E; aZawBU.:  
ez(i,1)=0; pe!dm}!h[  
ez(i,JE)=0; s7.p$r  
end; D-o7yc"K  
^0`<k
 
for i=1E; % 为了使每个磁场周围都有电场进行数组下标处理 EOBs}M
;  
for j=1:JE-1; /J[H5uA  
curl_e=ez(i,j)-ez(i,j+1); ,h@R' f!  
ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; <nb3~z1  
hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j)); &G pA1  
end; d@D;'2}Yc  
end; % 磁场HX循环结束 L

*UV  
`j}_BW_  
for i=1E-1; % 为了使每个磁场周围都有电场进行数组下标处理 aucZJjH  
for j=1:JE; }dd k}wga  
curl_e=ez(i+1,j)-ez(i,j); W) 33;E/}  
ihy(i,j)=ihy(i,j)+fj1(j)*curl_e; I NPYJ#%  
hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j)); oxMUW<gYd  
end; Pn+IJ=0Y  
end; % 磁场HY循环结束 G+hF [b44'  
end; V{T{0b"\U  
end;

谢谢大家,我找到问题了,系数弄错了.

怎么我看这个系数跟书上的一样呢？ 汗。。。我也遇到同样的问题。 as>:\hjP##  
还请lz说清楚下，问题在哪了？

学习了 pjN:&#Y]