三维pml matlab 程序，有问题，大侠帮忙啊，急！

超过48小时，复制至相关版块，此帖锁定
参与讨论

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。 9U~sRj=D  
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！ TGu]6NzyZ  
clear txXt<]N  
c=3e8;   b#E!wMClS  
muz=4*pi*1e-7;    0@dN$e  
epsz=1/(c*c*muz);       f3HleA&&  
f=1e+9;               !0 -[}vvU  
lambda=c/f;             mY,t]#^m7  
omega=2*pi*f; }?XNA.Wz  
ds=lambda/12;             X3,+aL`  
dt=ds/(2*c); ]~-
vU{  
muz=4*pi*1e-7;   "4e{Cq  
R0=1e-5;   )m[dfeqd +  
D'3. T{*rH  
%************************************************ R3Ka^l8R|  
ie=50;                                                                    %x向空间步             c-hhA%@Wq  
je=50;                                                                    %y向空间步       )90K^$93"  
ke=50;                                                                    %z向空间步 Ba /^CS  
ib=ie+1; m kHcGB!~  
jb=je+1; ;gNoiAxW  
kb=ke+1; j9/Ev]im|F  
is=25;                                                                     %激励源x向位置 C#4/~+  
js=25;                                                                     %激励源y向位置 RT
Ri{p  
ks=25;                                                                     %激励源z向位置 MJ0UZxnl  
iebc=8;                                                                    %左、右PML层层数 dt|f4XWF  
jebc=8;                                                                    %上、下PML层层数   \cRe,(?O  
kebc=8;                                                                    %前、后PML层层数
#Ez+1  
ibbc=iebc+1; _4#&!b6  
jbbc=jebc+1; g?d*cwtU  
kbbc=kebc+1; ;f)o_:(JJ  
nmax=100;                                                                  %时间步 bjYaJtn  
dlta=6;                                                                    %观测点距激励源步数 ZP^7`q)6  
>*cg K}!@  
*oU-V#  
%************************************************ B!zqvShF  
Exy=zeros(ie,jb,kb); 6%RN-  
Exz=zeros(ie,jb,kb); m5rJY/  
Eyx=zeros(ib,je,kb); @O%d2bgEWV  
Eyz=zeros(ib,je,kb); @%sr#YqY  
Ezx=zeros(ib,jb,ke); p7kH"j{xD  
Ezy=zeros(ib,jb,ke); hpOUz%  
Hxy=zeros(ib,je,ke); \w+a Q?e_  
Hxz=zeros(ib,je,ke); f:|O);nM  
Hyx=zeros(ie,jb,ke); FkJ
X)  
Hyz=zeros(ie,jb,ke); zS+_6s  
Hzx=zeros(ie,je,kb); Ar VNynQ  
Hzy=zeros(ie,je,kb); &i%1\o  
a#G]5TZ  
@ME .  
fid1=fopen('m1.txt','wt'); @CU~3Md*  
fid2=fopen('m2.txt','wt'); /l@h[}g+d-  
%************************************************ %:WM]dc  
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率 U?d4 ^  
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率 JDy;Jb  
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c; 0[T>UEI?  
!@vM@Z"  
nlkQ'XGAI  
for i=1:ib                                                                 %x向网格节点处电导率分布 s55t>t,g6  
                if      i<=iebc c/\$AJV.H  
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层 zr5(nAl  
                elseif i>=ib-iebc+1                                           %右侧PML层 xMAb=87_  
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2; {#'M3z=  
                else l`A4)8Y@  
                    sigx(i)=0;                                              %中间自由空间 `[z
d  
                end     q=40l  
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz)); K0Zq)<  
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz)); p3YF  
       Z2x%  
end   d$(>=gzBQ  
for i=1:ie                                                                 %x方向网格边中点处电导率分布 b[_${in:  
                if      i<=iebc Lc:DJA  
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层 SJdi*>  
                elseif i>=ie-iebc+1                                           %右侧PML层 eX@7f!uz  
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2; 2>bV+[@B  
                else @ dF]X  
                    sigxh(i)=0;                                              %中间自由空间 9M$N>[og  
                end                 /P3s.-sL  
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率       U\\nSU  
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz)); /K@{(=n  
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz)); 3!V$fl0  
end   p>@S61 & [  
9?_ybO~Oq  
6Y[|xu:N8Y  
wuQ>|\Zs  
for j=1:jb                                                                 %y方向网格节点处电导率分布 q4rDAQyPO  
                if      j<=jebc _p`@/[(|  
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层 2og8VI  
                elseif j>=jb-jebc+1                                           %上侧PML层 'o*:~n  
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2; e;/C}sK:  
                else G]- wN7G  
                    sigy(j)=0;                                              %中间自由空间 \1p5$0z  
                end               7berkU
0P  
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz)); 'F[ C 4  
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz)); 4fCg{  
end   L!]~J?)  
for j=1:je                                                                 %y方向网格边中点处电导率分布                 g}]EIv{  
                if      j<=jebc ttK`*Ng  
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层 '#b7Z?83C  
                elseif j>=je-jebc+1                                           %上侧PML层 Jqt&TqX@s  
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2; X;?Z_3I:5  
                else ,LHQ@/}A C  
                    sigyh(j)=0;                                              %中间自由空间 yw(E}   
                end 6Q6l?!|W4  
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率       V^U1o[`  
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz)); A|esVUo<3^  
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));         !&Z,ev  
end   BOQeP/>  
!dW77kLTg  
X0.-q%5  
H~P"uYKIZ  
for k=1:kb                                                                 %z方向网格节点处电导率分布 Fc"&lk4e  
                if      k<=kebc *KXg;777  
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层 F|DKp[<]8  
                elseif k>=kb-kebc+1                                           %后侧PML层 Twj?SV  
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2; #Rkldv'  
                else ;I+"MY7D  
                    sigz(k)=0;                                              %中间自由空间 6!3Jr  
                end oYG].PC  
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz)); P;dp>jL  
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));       9A4h?/  
end   l?/.uNw  
for k=1:ke                                                                 %z方向网格边中点处电导率分布                 >Lo!8Hen  
                if      k<=kebc `=0J:  
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层 XT|!XC!|  
                elseif k>=ke-kebc+1                                           %后侧PML层 r_kw "9  
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2; "k${5wk#Fl  
                else z]
YP  
                    sigzh(k)=0;                                              %中间自由空间 R;XR?59:.  
                end yT>t[t60/S  
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率     ^3-Wxn9&  
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz)); R "&(Ae?LR  
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));         DJ9;{,gm  
end   epH48)2  
/~LXY<-(  
7Pc0|Z/  
% Ex=Exy+Exz; ,ZYj8^gF  
% Ey=Eyx+Eyz; P9p{j1*;  
% Ez=Ezy+Ezx; g1uqsqYt  
| 3`qT#p{  
%************************************************ 47iwb  
tic Hhbf9)  
Qjj:r~l  
schedule_box=waitbar(0,'Please wait...'); % display the schedule box 7\ <4LX  
                                                                       Y"uFlHN&i  
for n=1:nmax     ,Dz2cR6  
.=>T yq  
           for i=1:ib $<UX/a\sH  
               for j=1:je UEq;}4Bo  
                 for k=1:ke e<>Lr  
                    if(i==is&&j==js&&k==ks) sX&M+
'h  
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 #&1Y!kbdd  
                    else !_`T8pJ`  
                        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)); id-VoHdK  
                    end Hr$oT=x[  
                                         wjmZ`UMz  
                 end Zw5\{Z0  
               end lK^Q#td:`  
           end /)i)wxi  
             sK:,c5^  
           T2wn!N?r  
            for i=1:ib 8i;N|:WdH  
               for j=1:je X*~NE\  
                 for k=1:ke n/"T7Y\2  
                   if(i==is&&j==js&&k==ks) '^l/e: (H3  
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 ]kmOX  
                   else 1q!JpC^  
                       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)); c=2e?  
                   end $p4aNC  
%                    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)); J{mP5<8>b  
               end ~^.&nph  
               end U
ZdE^Q[  
          end wS2iyrIB  
%           Hx=Hxy+Hxz; Y\T*8\h_[  
          for i=1:ie lO9{S=N  
              for j=1:jb x~GV#c  
                 for k=1:ke j8os6I  
                   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)); ru`;cXa,  
                 end =3=KoH/'  
               end Yi[dS`,d  
           end mm=Y(G[_%y  
          for i=1:ie                                                 s\&_Kbw]c  
               for j=1:jb LeW.uh3.  
                 for k=1:ke 4[3T%jA  
                   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)); Kfj*uzKB  
               end -aNTFt~|[  
               end ]tZ5XS  
          end 4(4JQ(5  
%           Hy=Hyx+Hyz; HHZ!mYr  
          for i=1:ie                                                   e>t9\vN#bx  
               for j=1:je EGwY|+3  
                 for k=1:kb ^c]Sl  
                   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)); ABV\:u  
               end vc2xAAQ  
               end yj$S?B Ee  
           end u,F d[[t  
          for i=1:ie                                                   1A-8,)  
               for j=1:je 4Uf+t?U9  
                 for k=1:kb wc%  
                   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)); 7o z(hO~  
               end 6\,^MI  
               end Za!c=(5  
          end #|XEBOmsQ  
% Hz=Hzx+Hzy; [6S"iNiyKT  
^@L[0Z`  
          for i=1:ie                                                   MatC2-aV1  
               for j=2:je  % PEC at j=1 & j=jb ^TWN_(-@  
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同
*vhm  
                   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)); "|DR"rr'j  
                 end !BEOeq@2.  
               end cM4?Ggn  
          end 16N8h]l  
          for i=1:ie                                               b=T+#Jb  
               for j=2:je =zA=D.D2  
                 for k=2:ke +$>ut r  
                   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)); FA^x|C=$  
                 end UKK}$B  
               end YLd 5  
          end ~v>w%]  
|/
fbU_d  
            for i=2:ie                                                     , m|9L{  
               for j=1:je ,:\2Lf  
                 for k=2:ke na']{a1K  
                   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)); 8 (^2  
                 end `A <
yDy  
               end r-YQsu&  
            end *gz{:}NX  
             [bPE?_a,  
         .. vB4cdW 2#3  
W,{`)NWg  

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

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

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

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

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

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

如果都没有迭代的话，自己调试一下吧。 BYr_Lz|T  
Hxy(is,js,ks)是源，那么Eyz和Ezy应该是可以得到值的，你把这几个分量都实时显示出来，看看问题到底在哪里吧。 $K6?(x_  
#!8^!}
nFO  
或者你给我个email，我发一个2D+PML的源码给你参考一下。

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

超过48小时，复制至相关版块，此帖锁定 8fJR{jD(s  
m.1LxM$8  
参与讨论