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

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。 #( 146  
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！ }Sh?S]]`  
clear mLLDE;7|}  
c=3e8;   l L@XM2"  
muz=4*pi*1e-7;   ,w:U#r~s"  
epsz=1/(c*c*muz);       eF-."1  
f=1e+9;               eiaF
aYe\  
lambda=c/f;             B!L{  
omega=2*pi*f; 9w"4K.  
ds=lambda/12;             "CQa.%  
dt=ds/(2*c); Vd+T$uC  
muz=4*pi*1e-7;   ['tY4$L(  
R0=1e-5;   Pw`8Wj  
e)? .r9pA;  
%************************************************ F8,RXlGfA[  
ie=50;                                                                    %x向空间步             }-2 2XYh  
je=50;                                                                    %y向空间步       `%"\@<  
ke=50;                                                                    %z向空间步 u[=r,^YQ  
ib=ie+1; xHLlMn4M  
jb=je+1; q\4Xs$APq  
kb=ke+1; bI9~jWgGp  
is=25;                                                                     %激励源x向位置 u.m[u)HQ  
js=25;                                                                     %激励源y向位置 A&Usddcp  
ks=25;                                                                     %激励源z向位置 oDAXiY$u  
iebc=8;                                                                    %左、右PML层层数 .2Elr(&*h  
jebc=8;                                                                    %上、下PML层层数   aP
@N)"  
kebc=8;                                                                    %前、后PML层层数 ?_9  
ibbc=iebc+1; LxSpctiNx  
jbbc=jebc+1; 2*l/3VW  
kbbc=kebc+1; 9ZsVy
 
nmax=100;                                                                  %时间步 T<Z &kYU:R  
dlta=6;                                                                    %观测点距激励源步数 0#Y5_i|p  
x}I+Iggi  
3J|F?M"N7  
%************************************************ :zke %Yx  
Exy=zeros(ie,jb,kb); `MN4uC  
Exz=zeros(ie,jb,kb); sfugY(m  
Eyx=zeros(ib,je,kb); By",rD- r  
Eyz=zeros(ib,je,kb); RmeD$>7  
Ezx=zeros(ib,jb,ke); A>
;bHf@  
Ezy=zeros(ib,jb,ke); Ed df2;-.  
Hxy=zeros(ib,je,ke); Z4w!p?Wqa  
Hxz=zeros(ib,je,ke); [ =9T*Sp  
Hyx=zeros(ie,jb,ke); j[G  
Hyz=zeros(ie,jb,ke); ;)z:fToh  
Hzx=zeros(ie,je,kb); Y0dEH^I  
Hzy=zeros(ie,je,kb); WH@,kH@  
U-(01-  
E`usknf>l  
fid1=fopen('m1.txt','wt'); h-K_Lr]  
fid2=fopen('m2.txt','wt'); m6\E$;`  
%************************************************ +RMSA^  
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率 -[9JJ/7y  
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率 Q}K"24`=  
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c; b;W3j  
&P}_bx  
}Gm>`cw-  
for i=1:ib                                                                 %x向网格节点处电导率分布 x$.^"l-vX  
                if      i<=iebc )9'K($  
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层 :tB1D@Cb6  
                elseif i>=ib-iebc+1                                           %右侧PML层 ;dtA4:IRZ4  
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2; l<LP&  
                else 03qQ'pq  
                    sigx(i)=0;                                              %中间自由空间 %i9E @EV  
                end     (
Ag16  
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz)); N06OvU2>xU  
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz)); y@:h4u"3  
       r,1!?s^L  
end   K6/Q}W
 
for i=1:ie                                                                 %x方向网格边中点处电导率分布 lH x^D;m6  
                if      i<=iebc FQ7T'G![  
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层 ):68%,  
                elseif i>=ie-iebc+1                                           %右侧PML层  vzs)[AD  
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2; 5z8d} I  
                else n&;85IF1  
                    sigxh(i)=0;                                              %中间自由空间 z2_*%S@  
                end                 HIR~"It$  
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率       vkx7paY_  
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz)); WwBOM~/`2  
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz)); #@9/g  
end   j@U]'5EVB  
F^t DL:  
]7F=u!/`<C  
L~rBAIdD  
for j=1:jb                                                                 %y方向网格节点处电导率分布  gmO!  
                if      j<=jebc Is)u }  
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层 gx8ouOh  
                elseif j>=jb-jebc+1                                           %上侧PML层 $%CF8\0  
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2; ]}-7_n#cC  
                else wOEj)fp.  
                    sigy(j)=0;                                              %中间自由空间 F|o:W75  
                end               :bu/^mW[  
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz)); 3G)#5Lf<  
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz)); 9Zt`u,;  
end   jrlVvzZ  
for j=1:je                                                                 %y方向网格边中点处电导率分布                 %S@ZXf~:  
                if      j<=jebc Pg0x/X{t  
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层 Jr ,;>   
                elseif j>=je-jebc+1                                           %上侧PML层 tqvN0vY5  
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2; a}B
Yov  
                else h-#6av:  
                    sigyh(j)=0;                                              %中间自由空间 7$vYo _  
                end hOu3 bA  
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率       ?qLFaFt/  
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz)); <ro7vPKNa  
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));         z0p*Z&  
end   iwZPpl";  
"3)C'WlEy/  
Pmr5S4Ka  
B:;pvW]  
for k=1:kb                                                                 %z方向网格节点处电导率分布 pMx*F@&nU  
                if      k<=kebc nJG U-Z  
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层 (
cAIvgI  
                elseif k>=kb-kebc+1                                           %后侧PML层 fcRxp{*zO  
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2; 3s,g*  
                else G_3O]BMKd)  
                    sigz(k)=0;                                              %中间自由空间 1#V_Z^OL  
                end zl>nSndRE  
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz)); hYT0l$Ng  
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));       av}k)ZT_  
end   fo*2:?K&  
for k=1:ke                                                                 %z方向网格边中点处电导率分布                 ]YnD  
                if      k<=kebc w;[NH/A^a  
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层 QuF:p  
                elseif k>=ke-kebc+1                                           %后侧PML层 w(*vj  
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2; 6y%qVx#!  
                else 7 S#J>*  
                    sigzh(k)=0;                                              %中间自由空间 L3u&/Tn2  
                end dUeN*Nq&(,  
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率     2\A$6N;_  
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz)); N ,'GN[s  
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));         %Q__!D[  
end   axv>6k  
N=T<_`$5  
RE7?KR>  
% Ex=Exy+Exz; }{K) 4M  
% Ey=Eyx+Eyz; W7R<%?  
% Ez=Ezy+Ezx; ;v)JnbsH}  
:t"^6xt  
%************************************************ ^e2VE_8L  
tic Xy|So|/bKd  
_wbF>z  
schedule_box=waitbar(0,'Please wait...'); % display the schedule box F'={q{2wH  
                                                                       jH5 k  
for n=1:nmax     Xk~D$~4<  
knu,"<  
           for i=1:ib oo/qb`-6  
               for j=1:je w=0(<s2  
                 for k=1:ke -j#2}[J7  
                    if(i==is&&j==js&&k==ks) &Cq`Y !y  
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 j\[dx^\=  
                    else iTBx\u%{  
                        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)); Uu10)/.LC  
                    end b$d;Qx  
                                         'Vzp2  
                 end 7{e 4c  
               end sQUM~HD\a  
           end ?Ny9'g>?  
             xBThq?N?  
           %B2'~|g  
            for i=1:ib fa jGZyd0:  
               for j=1:je E<{R.r  
                 for k=1:ke >a!/QMh  
                   if(i==is&&j==js&&k==ks) <$A  
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 h0*!;Z7  
                   else
aD<A.Lhy  
                       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)); Go`vfm"S  
                   end )Ys x}vSZ  
%                    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)); #px+;k5  
               end %~O,zs.2p  
               end ,8S/t+H  
          end ,uSMQS-O'4  
%           Hx=Hxy+Hxz; d\&U*=  
          for i=1:ie |k )=0mCz  
              for j=1:jb (Z+.45{-  
                 for k=1:ke s[RAHU  
                   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)); dc+>m,3$  
                 end }/0X'o  
               end 7X`g,b!  
           end <,3
a3  
          for i=1:ie                                                 1*P~!2h  
               for j=1:jb du $:jN\}  
                 for k=1:ke CmP9Q2  
                   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)); L4@K~8j7  
               end & TCkpS  
               end xeg/
A}yE  
          end )nC]5MXU  
%           Hy=Hyx+Hyz;
@+&LYy72  
          for i=1:ie                                                   GL>O4S<`  
               for j=1:je 9a[9i}_  
                 for k=1:kb >2y':fO  
                   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)); 5N#aXG^9  
               end Q
GMV}y  
               end G*?8MTP8![  
           end ~dyTVJ$  
          for i=1:ie                                                   oM
X  
               for j=1:je b<tNk]7  
                 for k=1:kb fF!Yp iI"  
                   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)); ]{;gw<T  
               end $g^@AdE%  
               end wm+};L&_  
          end _>o:R$ %}  
% Hz=Hzx+Hzy; HDz5&7* .  
z 4e7PW|  
          for i=1:ie                                                   AmUr.ofu  
               for j=2:je  % PEC at j=1 & j=jb SpIv#?  
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同 U45e2~1!O  
                   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)); lB8-Z ow  
                 end @f_Lp%K  
               end BwN0!lsF3  
          end ")1:F>  
          for i=1:ie                                               XnH05LQ  
               for j=2:je o3XvRj  
                 for k=2:ke ^)470K`%)  
                   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)); :p1u(hflS  
                 end H9Gh>u]}  
               end m;$b'pT  
          end pN,u
`[  
^Y?k0z  
            for i=2:ie                                                     G~]Uk*M q  
               for j=1:je k`cfG\;r  
                 for k=2:ke <X#C)-.  
                   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)); K6)Gc%:`  
                 end vRTkgH#4l  
               end }g@v`5  
            end c^5~QGuQ  
             \:P>l
e'1  
           .. P7ao5NP  
CmWeY$Jb  

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

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

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

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

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

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

如果都没有迭代的话，自己调试一下吧。 ZECfR>`x  
Hxy(is,js,ks)是源，那么Eyz和Ezy应该是可以得到值的，你把这几个分量都实时显示出来，看看问题到底在哪里吧。 fJg+Ryo  

H:|uw  
或者你给我个email，我发一个2D+PML的源码给你参考一下。

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

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

你是直接做的3D的吗？