[已解决]三维pml matlab 程序，有问题，大侠帮忙啊，急！

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。 zOiu5  
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！ (UxW;  
clear V=*wKuB  
c=3e8;   &oX>*6L  
muz=4*pi*1e-7;   RVQh2'w  
epsz=1/(c*c*muz);       X)% A6M  
f=1e+9;               WILMH`  
lambda=c/f;             q?8| [.  
omega=2*pi*f; % S os  
ds=lambda/12;             {Ja!~N;3  
dt=ds/(2*c); v'3J.?N  
muz=4*pi*1e-7;   Dbz3;t  
R0=1e-5;   zl
d#qG6  
a5TioQ  
%************************************************ q0zr E5  
ie=50;                                                                    %x向空间步             sjV!5Z  
je=50;                                                                    %y向空间步       gp\<p-}  
ke=50;                                                                    %z向空间步 BGX.U\uc  
ib=ie+1; K9up:.{QQ  
jb=je+1; Kuu *&u  
kb=ke+1; uwy:t!(j  
is=25;                                                                     %激励源x向位置 M"94#.dKK  
js=25;                                                                     %激励源y向位置 RtM8yar+sn  
ks=25;                                                                     %激励源z向位置 w {3<{  
iebc=8;                                                                    %左、右PML层层数 =aTv! 8</  
jebc=8;                                                                    %上、下PML层层数   *vwbgJG! *  
kebc=8;                                                                    %前、后PML层层数 av|g}xnj  
ibbc=iebc+1; 1bn^.768l  
jbbc=jebc+1; FFEfI4&SfS  
kbbc=kebc+1; Zo~  
nmax=100;                                                                  %时间步 \r+8qC[,  
dlta=6;                                                                    %观测点距激励源步数 e0,|Wm  
lp^<3o*1  
fd.^h*'mU  
%************************************************ pWJFz-  
Exy=zeros(ie,jb,kb); "L&#lfOKG  
Exz=zeros(ie,jb,kb); %Da1(bBh  
Eyx=zeros(ib,je,kb); 0a XPPnuX  
Eyz=zeros(ib,je,kb); z+n,uHs  
Ezx=zeros(ib,jb,ke); VG ;kPzze  
Ezy=zeros(ib,jb,ke); ~G6Ox)/  
Hxy=zeros(ib,je,ke); bl&nhI)w  
Hxz=zeros(ib,je,ke); ][p>Y>:b-  
Hyx=zeros(ie,jb,ke); z\%67C  
Hyz=zeros(ie,jb,ke); m3/O.DY%0  
Hzx=zeros(ie,je,kb); )`O~f_pIC  
Hzy=zeros(ie,je,kb); #6HA\dE  
7V!*NBsl  
|$ lM#Ua  
fid1=fopen('m1.txt','wt'); Rx=>6,)'  
fid2=fopen('m2.txt','wt'); o9dY9o+Z  
%************************************************ 4?q
<e*W  
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率 ?b>,9A.Z  
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率 2OVRf0.R~  
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c;
_v> }_S  
EZ`te0[  
fy@<&U5rg  
for i=1:ib                                                                 %x向网格节点处电导率分布 6"&6`f  
                if      i<=iebc S(*sw 0O@+  
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层 hFy;ffs.  
                elseif i>=ib-iebc+1                                           %右侧PML层 ?q{,R"  
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2; %UERc{~o*,  
                else eEv@}1~  
                    sigx(i)=0;                                              %中间自由空间 :Ra,Eu  
                end     ?WqT[MnK  
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz)); `3WFjU5a  
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz)); (Hb:?(  
       +{f:cea (1  
end   lHPd"3HDK  
for i=1:ie                                                                 %x方向网格边中点处电导率分布 n|RJ;d30Q  
                if      i<=iebc q%"VYt4  
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层 CU@Rob}s  
                elseif i>=ie-iebc+1                                           %右侧PML层 '9 [vDG~  
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2; 9CWezI+  
                else 9\mLW"  
                    sigxh(i)=0;                                              %中间自由空间 _n50C"X=&(  
                end                 ?En O"T.  
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率       zGkS^Z=(  
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz)); Gsq00j &<Z  
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz)); 8h*Icf  
end   '
O_3)x5  
m4hg'<<V  
}o?APvd  
SVh 7zh  
for j=1:jb                                                                 %y方向网格节点处电导率分布 ($;77fPR  
                if      j<=jebc O @j} K4  
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层 &-Gqdnc  
                elseif j>=jb-jebc+1                                           %上侧PML层 PIoLywpRn  
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2; eo?;`7  
                else _4U5  
                    sigy(j)=0;                                              %中间自由空间 %u5L!W&  
                end               Kzm+GW3o[  
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz)); 2j}\3Pi  
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz)); xRzFlay8  
end   e]$}-i@#  
for j=1:je                                                                 %y方向网格边中点处电导率分布                 (mTE;s(  
                if      j<=jebc YA_c N5p/@  
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层 `tA" }1;ka  
                elseif j>=je-jebc+1                                           %上侧PML层 g+
Sbl  
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2; qzKdQ&vO  
                else g4=pnK8  
                    sigyh(j)=0;                                              %中间自由空间 i%#+\F.&  
                end *U,@q4  
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率       ;S^
'V  
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz)); 2a`o &S  
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));         SwTL|+u  
end   "!ug_'VW  
d"*uB
VzXm  
v4`"1Ss,K  
#FCnA  
for k=1:kb                                                                 %z方向网格节点处电导率分布 Mb|a+,:>3  
                if      k<=kebc |@ s,XS  
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层 ;5S9y7[i|  
                elseif k>=kb-kebc+1                                           %后侧PML层 :OhHb#D  
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2; t\k$};qJ  
                else g9}Dn
CT*.  
                    sigz(k)=0;                                              %中间自由空间 -w}]fb2Q>  
                end p6#g;$V$  
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz)); Eg#K.5hJ  
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));       _/-jX  
end   z<U-#k7nz  
for k=1:ke                                                                 %z方向网格边中点处电导率分布                 Oj3.q#)`Z  
                if      k<=kebc /v1Q4mq  
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层 P3x= 8_#  
                elseif k>=ke-kebc+1                                           %后侧PML层 =S+wCN  
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2; _VRpI)mu  
                else mD$A4Y-'p  
                    sigzh(k)=0;                                              %中间自由空间 1{ ~#H<K  
                end N~goI#4  
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率     hKLCJ#
T  
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz)); }Qn&^[[miL  
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));         @"Fme-~
 
end   mS$j?>m  
cdl&9-}  
.f%fHj  
% Ex=Exy+Exz; *`ua'"="k  
% Ey=Eyx+Eyz; n22zq6m  
% Ez=Ezy+Ezx; =jOv] /  
"U>JM@0DNm  
%************************************************ t{^*6XOcJ  
tic (2J: #  
6}[I2F_^  
schedule_box=waitbar(0,'Please wait...'); % display the schedule box [!HEQ8 2g  
                                                                       5\5/   
for n=1:nmax     PV'x+bN5  
=.f-w0V  
           for i=1:ib r@h5w_9  
               for j=1:je #~}nFY.  
                 for k=1:ke Wuc S:8#|  
                    if(i==is&&j==js&&k==ks) ZM!CaR  
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 oTU!R ,  
                    else jnKWZ/R  
                        Hxy(i,j,k)=cpy(j).*Hxy(i,j,k)+cqy(j).*(Ezx(i,j,k)+Ezy(i,j,k)-Ezx(i,j+1,k)-Ezy(i,j+1,k)); |[<_GQl  
                    end U@_dm/;0&  
                                         \o}xF@sM5  
                 end z;{iM/Xe  
               end TN!j13,  
           end U\4g#!qj  
             `#F{Waww'  
           g]<4&)~  
            for i=1:ib d6}r#\  
               for j=1:je D0&,?  
                 for k=1:ke Z0x ar]4V  
                   if(i==is&&j==js&&k==ks) Z&Pg"a?\  
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 bH7X'%r  
                   else jVv0ST*z  
                       Hxz(i,j,k)=cpz(k).*Hxz(i,j,k)+cqz(k).*(Eyx(i,j,k+1)+Eyz(i,j,k+1)-Eyx(i,j,k)-Eyz(i,j,k)); bf ]f=;.+  
                   end A-Sv;/yD_  
%                    Hxz(i,j,k)=cpz(k).*Hxz(i,j,k)+cqz(k).*(Eyx(i,j,k+1)+Eyz(i,j,k+1)-Eyx(i,j,k)-Eyz(i,j,k)); $2oTkOA   
               end bhTb[r  
               end D.B.7-_8  
          end &zl|87M  
%           Hx=Hxy+Hxz; R} eN@#"D  
          for i=1:ie :q$.,EZ4#n  
              for j=1:jb >Ea8G,  
                 for k=1:ke
-BrMp%C  
                   Hyz(i,j,k)=cpz(k).*Hyz(i,j,k)+cqz(k).*(Exy(i,j,k)+Exz(i,j,k)-Exy(i,j,k+1)-Exz(i,j,k+1)); nhB
1D-  
                 end YSr9VpqWV  
               end Y;dz,}re  
           end 0bceI  
          for i=1:ie                                                 GY6`JWk  
               for j=1:jb \\PjKAsh  
                 for k=1:ke f=(?JT  
                   Hyx(i,j,k)=cpx(i).*Hyx(i,j,k)+cqx(i).*(Ezx(i+1,j,k)+Ezy(i+1,j,k)-Ezx(i,j,k)-Ezy(i,j,k)); }iXDa?6%  
               end |%F=po>w  
               end W98i[Q9A7  
          end s:>VaGC  
%           Hy=Hyx+Hyz; "Gfh,e  
          for i=1:ie                                                   bR*-Ht+wd  
               for j=1:je l4D+Y  
                 for k=1:kb {@H6HqD  
                   Hzx(i,j,k)=cpx(i).*Hzx(i,j,k)+cqx(i).*(Eyx(i,j,k)+Eyz(i,j,k)-Eyx(i+1,j,k)-Eyz(i+1,j,k)); jqWu  
               end g`{;(/M+  
               end ^crCy-`#  
           end x5,++7Tz  
          for i=1:ie                                                   C]O(T2l{l  
               for j=1:je lGV0*Cji  
                 for k=1:kb rHC>z7+z.  
                   Hzy(i,j,k)=cpy(j).*Hzy(i,j,k)+cqy(j).*(Exy(i,j+1,k)+Exz(i,j+1,k)-Exy(i,j,k)-Exz(i,j,k)); Q3n,)M[N  
               end %+@O#P  
               end Yl4^AR&  
          end m'Amli@[  
% Hz=Hzx+Hzy; 9TgIB  
~bM4[*Q7  
          for i=1:ie                                                   _GXk0Ia3`  
               for j=2:je  % PEC at j=1 & j=jb N=4G=0 `ke  
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同 b*;Si7-  
                   Exy(i,j,k)=cay(j).*Exy(i,j,k)+cby(j).*(Hzx(i,j,k)+Hzy(i,j,k)-Hzx(i,j-1,k)-Hzy(i,j-1,k)); ^1S!F-H4\  
                 end v~f HYa>  
               end V, Z|tB^  
          end <{dVKf,e  
          for i=1:ie                                               0[RL>;
D:  
               for j=2:je h;C5hU4P  
                 for k=2:ke [;r)9mh7  
                   Exz(i,j,k)=caz(k).*Exz(i,j,k)+cbz(k).*(Hyx(i,j,k-1)-Hyx(i,j,k)+Hyz(i,j,k-1)-Hyz(i,j,k)); Ttu2skcv  
                 end Sz%tJD..  
               end G"-?&)M#a  
          end yQ_B)b  
mC4zactv  
            for i=2:ie                                                     &vo--V1|  
               for j=1:je @oNH@a j%  
                 for k=2:ke Buf/@B7+\  
                   Eyz(i,j,k)=caz(k).*Eyz(i,j,k)+cbz(k).*(Hxy(i,j,k)-Hxy(i,j,k-1)+Hxz(i,j,k)-Hxz(i,j,k-1));  ,V,`Jf  
                 end :8L8q<U  
               end Jv>gwV{  
            end bx#>BK!  
             PXK7b2fE.  
           .. ;;_,~pI?k  
i2@VB6]?  

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

你说电磁场迭代有大问题，到底是什么大问题，是电磁波传播不稳定？很快就变成无穷大了么？还是计算结果和预期不对？

是除了我的激励那有波形，其他点一律没有产生迭代

第一是PML煤质中损耗电导率的计算在不同的场量(相差半网格)上是不一样的。

分离场的PML，是最繁琐的、代码量最多的PML方法，大体一看代码的长度，就知道算法编写有问题。

分离的场步进计算的系数看着有点别扭啊

如果都没有迭代的话，自己调试一下吧。 ZO4*sIw%  
Hxy(is,js,ks)是源，那么Eyz和Ezy应该是可以得到值的，你把这几个分量都实时显示出来，看看问题到底在哪里吧。 GN!qyT  
$BFvF ,n  
或者你给我个email，我发一个2D+PML的源码给你参考一下。

碟代/碟代/碟代/碟代/

版主,最近一直没上网,不好意思,这个程序我还在研究,下个帖子我发最新得程序,但仍未有最好效果,也就是说这个问题还没有解决掉.

你是直接做的3D的吗？