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

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。 x%> e)L<  
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！ mj,fp2D;%  
clear Wsj=!Obc  
c=3e8;   `TKD<&oL  
muz=4*pi*1e-7;    -7]Xjb5  
epsz=1/(c*c*muz);       T=Z.U$  
f=1e+9;               W 0[N0c  
lambda=c/f;             Uu p(6`7  
omega=2*pi*f; [bQj,PZ&  
ds=lambda/12;             b3qc_  
dt=ds/(2*c); ]gb=  
muz=4*pi*1e-7;   ^s/  
R0=1e-5;   9l/EjF^  
ZXb0Y2AVx  
%************************************************ g^>#^rLU  
ie=50;                                                                    %x向空间步             2(YPz|~W  
je=50;                                                                    %y向空间步       |FR3w0o  
ke=50;                                                                    %z向空间步 mryN}  
ib=ie+1; xycH~ ?  
jb=je+1; &~sfYW  
kb=ke+1; |BU+:+
 
is=25;                                                                     %激励源x向位置 j8,n7!G  
js=25;                                                                     %激励源y向位置 CZ{k@z
`r  
ks=25;                                                                     %激励源z向位置 tJI,r_  
iebc=8;                                                                    %左、右PML层层数 ?:rx1}:F  
jebc=8;                                                                    %上、下PML层层数   jdqj=Yc  
kebc=8;                                                                    %前、后PML层层数 /'DwfX  
ibbc=iebc+1; 4jmK].  
jbbc=jebc+1; *<V^2z$y_  
kbbc=kebc+1; )|#%Czd4  
nmax=100;                                                                  %时间步 |01?w|  
dlta=6;                                                                    %观测点距激励源步数 rJfqA@  
=v6*|  
``Q2P%  
%************************************************ ZFh+x@  
Exy=zeros(ie,jb,kb); b B  
Exz=zeros(ie,jb,kb); 2YwVU.*>  
Eyx=zeros(ib,je,kb); %,*$D}H  
Eyz=zeros(ib,je,kb); cd@.zg'sYn  
Ezx=zeros(ib,jb,ke); /i|z.nNO  
Ezy=zeros(ib,jb,ke); :K.4n  
Hxy=zeros(ib,je,ke); I=DLPgzO9  
Hxz=zeros(ib,je,ke); B)SLG]72f  
Hyx=zeros(ie,jb,ke); dx_6X!=.J  
Hyz=zeros(ie,jb,ke); 3eIr{xs  
Hzx=zeros(ie,je,kb); l0]
d  
Hzy=zeros(ie,je,kb); @!tVr3;N$  
-C7FuD[Xw  
;^k7zNf-  
fid1=fopen('m1.txt','wt'); S4NL "m  
fid2=fopen('m2.txt','wt'); LX+5|u  
%************************************************ JQKdW  
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率 B\`Aojw"E?  
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率 %Q|eiXD  
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c; /!FWuRe^  
7LZ^QC  
Al"3 kRJJ  
for i=1:ib                                                                 %x向网格节点处电导率分布 p QE)p  
                if      i<=iebc y7u^zH6wj  
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层 v@ C,RP9  
                elseif i>=ib-iebc+1                                           %右侧PML层 8|1^|B(l  
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2; i^Ut015q%  
                else YrI|gz)  
                    sigx(i)=0;                                              %中间自由空间 mPZGA\  
                end     ]5CNk+`'  
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz)); .q`{Dgc~  
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz)); `i!wq&1g7  
       V-O(U*]  
end   EYS>0Y  
for i=1:ie                                                                 %x方向网格边中点处电导率分布 jB$SUO`*  
                if      i<=iebc ;I6s-moq_  
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层 0`#(Toe{B  
                elseif i>=ie-iebc+1                                           %右侧PML层 |lQ;ALH!  
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2; #~ v4caNx  
                else `vk0c  
                    sigxh(i)=0;                                              %中间自由空间 2i=H"('G)+  
                end                 eBW=^B"y+  
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率       h#YD~!aJ  
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz)); ?^!,vh  
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz)); X!,@j\L  
end   8}p8r|d!ls  
aC=['a>)  
9>*c_  
Z)< wv&K  
for j=1:jb                                                                 %y方向网格节点处电导率分布 $r.U  
                if      j<=jebc Ibr%d2yS=  
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层 Xh3;  
                elseif j>=jb-jebc+1                                           %上侧PML层 4n( E;!s  
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2; F}?<v8#z0  
                else [W,|kDK  
                    sigy(j)=0;                                              %中间自由空间 !AN^ ,v]D  
                end               jRZ%}KX  
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz)); dJ3IUe  
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz)); N_wp{4 0/  
end   hXL|22>w<  
for j=1:je                                                                 %y方向网格边中点处电导率分布                 _ *(bmJM  
                if      j<=jebc 6.Ie\5-a;  
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层 %n?vJ#aX%  
                elseif j>=je-jebc+1                                           %上侧PML层 ;6;H*Y0,|E  
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2; h^|5|l  
                else {+T/GBF-K=  
                    sigyh(j)=0;                                              %中间自由空间 V-W'RunnW  
                end R'K/t|MC  
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率       3x*z\VJ  
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz)); m9$:9yRm  
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));         bi[7!VQf  
end   G$WOzY(  
l[h'6+o  
.-I|DVHe  
T|4snU2M  
for k=1:kb                                                                 %z方向网格节点处电导率分布 Zc5 :]]  
                if      k<=kebc `^F: -  
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层  !rL<5L  
                elseif k>=kb-kebc+1                                           %后侧PML层 'Z!Ga.I  
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2;
1i|.h  
                else Mhpdaos  
                    sigz(k)=0;                                              %中间自由空间 ^QL 877  
                end t%0?N<9YkU  
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz)); *YV S|6bs  
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));       o \L!(hm  
end   >l1r,/\\
 
for k=1:ke                                                                 %z方向网格边中点处电导率分布                 NlPS#  
                if      k<=kebc 0&w0aP`Y  
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层 v|t{1[C  
                elseif k>=ke-kebc+1                                           %后侧PML层 7} O;FX+x  
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2; M?&zY "c  
                else laREjN/\`  
                    sigzh(k)=0;                                              %中间自由空间 HMQI&Lh=U  
                end Opg#*w%-  
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率     J1F{v)T'?  
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz)); gT52G?-  
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));         m2bDHQ+  
end   u=K2Q4  
I44s(G1jl  
j~av\SCU*  
% Ex=Exy+Exz; Q
JX/7RA  
% Ey=Eyx+Eyz; a+z2Zd!u\x  
% Ez=Ezy+Ezx; tai Vk4  
7XAvd-  
%************************************************ ? Nj)6_&  
tic @PkJY  
b)}+>Wx  
schedule_box=waitbar(0,'Please wait...'); % display the schedule box C3;[e0.1b  
                                                                       ~1 ZD[@  
for n=1:nmax     pA4 ,@O  
8n)WW$  
           for i=1:ib @'jfKW  
               for j=1:je wT_^'i*@I  
                 for k=1:ke W4(?HTWZ  
                    if(i==is&&j==js&&k==ks) /'&;Q7!)  
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 yY=<'{!  
                    else G?)vWM`j  
                        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)); "ED8z|]j  
                    end ;Rnhe_A.  
                                         _n O.-  
                 end s@@1 *VQ  
               end bg7n  
           end "Ai6<:ml  
             ]|oJ)5P  
           DZ`k[Z.VZ  
            for i=1:ib ;)/
@Xx  
               for j=1:je D(XqyN-P  
                 for k=1:ke LJuW${Y  
                   if(i==is&&j==js&&k==ks) d0UZ+ RR#  
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 g`Md80*Zfk  
                   else KCqqJ}G  
                       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)); }[=YU%[o:  
                   end &"S/Lt  
%                    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));  W'/>et  
               end ;bjnL>eW  
               end \HSicV#i  
          end
9%,;XQ  
%           Hx=Hxy+Hxz; F?]N8W  
          for i=1:ie y?:dE.5p|  
              for j=1:jb YMzBAf  
                 for k=1:ke X>,A  
                   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)); 3)E(RyQA3  
                 end 78
Y@OL_$  
               end zJl_ t0  
           end 9YKDguG  
          for i=1:ie                                                 y:so L:(F  
               for j=1:jb JeQ[qQ  
                 for k=1:ke h}Ygb-uZ  
                   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)); 9D_wG\g  
               end d<T%`:s<  
               end $G";2(-k  
          end ~<O,Vs_C/  
%           Hy=Hyx+Hyz; :t$A8+A+0  
          for i=1:ie                                                   v*TeTA %  
               for j=1:je 3E|;r _; 8  
                 for k=1:kb {*xBm#  
                   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)); {BOLPE-  
               end ~ =.CTm]vf  
               end 2P}I'4C-  
           end FmT `Oa>  
          for i=1:ie                                                   :Yvb
U Y  
               for j=1:je 8&(-8  
                 for k=1:kb I'_v{k5ZI  
                   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)); aqL<v94wX  
               end #x)}29%e#  
               end L#j/0IHD  
          end f%Ke8'&  
% Hz=Hzx+Hzy; CT"Fk'B'  
d6b.z
P  
          for i=1:ie                                                   TmP8q  
               for j=2:je  % PEC at j=1 & j=jb TVkcDS  
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同 zO.6WJ  
                   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)); #93}E Y  
                 end j#igu#MB
*  
               end K2|7%  
          end V5s&hZZYa  
          for i=1:ie                                               *{[d%B<lp  
               for j=2:je FdxsUDL  
                 for k=2:ke DrS?=C@  
                   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)); I'A:J  
                 end #v/ry)2Y=  
               end yYX :
huw  
          end M,..Kw/ }~  
l2/@<0P  
            for i=2:ie                                                     fv?45f  
               for j=1:je {j@+h%sF>+  
                 for k=2:ke otnV-7)@  
                   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)); M&e8zS  
                 end VkXn8J  
               end }tR'Hz2  
            end #s^s_8#&e  
             -1mvhR~  
           .. .<w)Bmh  
Wem?{kx0  

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

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

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

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

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

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

如果都没有迭代的话，自己调试一下吧。 +N@F,3yNa  
Hxy(is,js,ks)是源，那么Eyz和Ezy应该是可以得到值的，你把这几个分量都实时显示出来，看看问题到底在哪里吧。 mhuaXbr  
[1GEe
 
或者你给我个email，我发一个2D+PML的源码给你参考一下。

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

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

你是直接做的3D的吗？